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The supermassive black holes (SMBHs) massive enough (> fewxlO^MQ) to 
^ . power the bright redshift z ^ 6 quasars observed in the Sloan Digital Sky Sur- 



00 



vey (SDSS) are thought to have assembled by mergers and/or gas accretion 

^ ! from less massive "seed" BHs. If the seeds are the ~ IO^Mq remnant BHs 

I— i! from the first generation of stars, they must be in place well before redshift 

p ,! z = 6, and must avoid being ejected from their parent proto-galaxies by the 

Q ! large (severalxlO^ km s~^) kicks they suffer from gravitational-radiation induced 

^ ! recoil during mergers with other BHs. We simulate the SMBH mass function at 

c^ , redshift z > 6 using dark matter (DM) halo merger trees, coupled with a pre- 

! scription for the halo occupation fraction, accretion histories, and radial recoil 

^ ! trajectories of the growing BHs. Our purpose is (i) to map out plausible scenar- 

CN ! ios for successful assembly of the z ^ 6 quasar BHs by exploring a wide region 

f~^ i of parameter space, and (ii) to predict the rate of low-frequency gravitational 

"^ [ wave events detectable by the Laser Interferometer Space Antenna (LISA) for 

^ I each such scenario. Our main findings are as follows: (1) ~ IOOMq seed BHs 

OO ! can grow into the SDSS quasar BHs without super-Eddington accretion, but 

.. I only if they form in minihalos at z > 30 and subsequently accrete > 60% of the 

.!^ I time; (2) the scenarios with optimistic assumptions required to explain the SDSS 

rN I quasar BHs overproduce the mass density in lower-mass (fewxlO^Mo < Mbh ^ 

5^ \ fewx IO^Mq) BHs by a factor of 10^ — 10^, unless seeds stop forming, or accrete at 

a severely diminished rates or duty cycles (e.g. due to feedback), at ^ ^ 20 — 30. 

We also present several successful assembly models and their LISA detection 

rates, including a "maximal" model that gives the highest rate (~ 30 yr~^ at 

z = 6) without overproducing the total SMBH density. 



Subject headings: cosmology: theory - galaxies: formation - quasars: general - 
black hole physics - accretion 
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1. Introduction 

The discovery of briglit quasars at redsliifts z > 6 in the Sloan Digital Sky Survey 
(SDSS) implies that black holes (BHs) as massive as several x IO^Mq were already assembled 
when the age of the universe was less than ^ 1 Gyr (see the recent review by Fan 2006). 
These objects are among the oldest detected discrete sources of radiation in the Universe. 
The likelihood that all of these quasars are significantly magnified by gravitational lensing, 
without producing detectable multiple images (Richards et al. 2004), is exceedingly small 
(Keeton et al. 2005), and if their luminosities are powered by accretion at or below the 
Eddington rate, the central objects must be ~ 10^ Mq supermassive black holes (SMBHs). 
In particular, the quasar SDSS J1148+5251 (Fan et al. 2003, 2001) is likely to be powered 
by a SMBH with a mass of ^ IO^-^Mq (Willott et al. 2003). 

The mechanism by which such massive BHs formed within 1 Gyr after the Big Bang 
remains poorly understood. Generically, these SMBHs are thought to have assembled by 
mergers with other BHs and/or by gas accretion [J onto less massive BHs. If the first ("seed") 
BHs are the ~ IO^Mq remnant BHs of the first generation of stars (e.g. Heger et al. 2003), 
they must be in place well before redshift z = Q. If accretion onto BHs is limited at the 
Eddington rate with radiative efficiency e, defined as the fraction of the rest mass energy of 
matter falling onto the BH that is released as radiation, then 1 — e of the matter is accreted 
and the growth of the BH mass m is given by 

dlnm 1 — eATcGfinip 
at e a^c 

where G is the gravitational constant, c is the speed of light, /i ~ 1.15 is the mean atomic 
weight per electron for a primordial gas, and Ue is the Thompson electron cross section. The 
e-folding time scale for mass growth is tEdd ~ 4.4 x 10^ yr for e = 0.1. In the concordance 
cosmological model (see below) the time elapsed between redshifts 2; = 30 (when the first 
seeds may form) and z = 6.4 (the redshift of the most distant quasar) is ~ 0.77 Gyr, allowing 
for a mass growth by a factor of ~ 10^'^. Therefore, individual ~ IOOMq seeds can grow into 
the SDSS quasar BHs through gas accretion alone, provided the accretion is uninterrupted 
at close to the Eddington rate and e ^ 0.1. A higher efficiency and/or a lower time-averaged 
accretion rate will require many seed BHs to merge together; the number of required mergers 
increases exponentially for lower time-averaged accretion rates. 

The discovery of the bright quasars at z > 6 were followed by the first successful 
numerical calculations in full general relativity of the coalescence of a BH binary and the 



^In this paper, "accretion" onto BHs should be assumed to mean gas accretion, unless otherwise noted. 
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corresponding emission of gravitational waves (GWs) (Pretorius 2005; Campanelli et al. 
2006; Baker et al. 2006). These calculations also confirmed a result previously known from 
post-Newtonian (Kidder 1995) and perturbation-theory treatments (Favata et al. 2004, 
Schnittman & Buonanno 2007): the coalesced product receives a large center-of-mass recoil 
imparted by the net linear momentum accumulated by the asymmetric gravitational wave 
emission (see Schnittman et al. 2008 for a recent detailed discussion of the physics of the 
recoil, and for further references). Typical velocities for this gravitational recoil (or "kick") 
are in excess of ~ 100 km s~^. This is likely more than sufficient to eject the BHs residing 
in the low-mass protogalaxies in the early Universe, since the escape velocities from the DM 
halos of these galaxies are only a few km s~^. 

Several recent works have studied the role of gravitational kicks as an impediment to 
merger-driven modes of SMBH assembly. Simple semi-analytic models show that if every 
merger resulted in a kick large enough to remove the seed BHs from halos with velocity 
dispersions up to ~ 50 km s~^, then super-Eddington accretion would be required to build 
SMBHs of the required mass in the available time (Haiman 2004; Shapiro 2005). Monte- 
Carlo merger tree models that exclude kicks entirely (Bromley et al. 2004) or which include 
a distribution of kick velocities extending to low values (e.g. Yoo & Miralda-Escude 2004 
[hereafter YM04]; Volonteri & Rees 2006 [hereafter VR06]) give a slightly more optimistic 
picture, showing that if seed BHs form in most minihalos in the early Universe, and especially 
if ejected seeds are rapidly replaced by new seeds (YM04), then SMBH assembly is just 
possible before z ^ 6 without exceeding the Eddington accretion rate. These works are 
encouraging steps toward demonstrating that there are plausible physical models that lead 
to the timely assembly of SMBHs massive enough to power the z > 6 SDSS quasars. 

At present, we have no direct observational constraints on SMBH assembly at 2; > 6, 
and there is, in principle, a large range of "physically plausible" possibilities. The Laser 
Interferometer Space Antenna (LISA) is expected to be able to detect mergers of SMBHs in 
the mass range ~ (10^-10^) M0/(1 + z) out to z ~ 30, and to extract binary spins and BH 
masses with high precision up to 2; ~ 10 (Vecchio 2004; Lang & Hughes 2006, 2007). It is also 
likely that by the time LISA is operational, there will be additional independent constraints 
on the demography of high-redshift SMBHs. It is therefore a useful exercise to calculate 
the expected LISA event rate from high-redshift SMBH mergers (e.g. Wyithe & Loeb 2003; 
Sesana et al. 2004, 2005, 2007) in a range of plausible models. Note that published estimates 
(Menou et al. 2001, Heger et al. 2003, Menou 2003, Haehnelt 2003,Wyithe & Loeb 2003, 
Sesana et al. 2004, Islam et al. 2004, Koushiappas & Zentner 2006, Micic et al. 2007, Lippai 
et al. 2008, Arun et al. 2008) for the LISA event rate, even at lower redshifts, vary by 
orders of magnitude, from ~ 1 to as high as ~ 10^ yr~^; there is a large raiige even among 
models that explicitly fit the evolution of the quasar luminosity function (I12!). A related 
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open question is to what degree the LISA data stream can help pinpoint the actual SMBH 
assembly scenario. One aim of this paper is to understand the model degeneracies that can 
lead to similar LISA data streams. Another is to explore as much as possible the full variety 
of LISA event rates arising from various "physically plausible" assembly models. 

The physical factors that determine the growth of SMBHs at high redshift fall broadly 
into four categories: (1) the properties of the initial seed BHs, such as their redshift, mass, 
and abundance; (2) the time-averaged gas accretion rate of individual seeds; (3) the merger 
rates of BHs; and (4) effects governing the fate of gravitationally kicked BHs. The first 
category determines the number and mass of seed BHs available for assembly, and depends 
primarily on the behavior of gas in the host DM halos, and the mass and metallicity of the first 
stars. The second category measures the subsequent growth through accretion, and depends 
on the availability of fuel over the Hubble time, and its ability to shed angular momentum 
and accrete onto the BH. The third category is a combination of the halo merger rate, the 
rarity of seeds, and the timescale for the formation, orbital decay, and ultimate coalescence 
of a SMBH binary. Finally, the recoil velocity of the coalesced binary is determined by 
the mass ratio and spin vectors of the BHs, and its subsequent orbit - and whether it is 
retained or ejected before the next merger - will be determined by the overall depth of the 
gravitational potential of the DM host halo, and on the spatial distribution of gas and DM 
in the central region of the halo, which determine drag forces on the kicked BH. 

Our approach to model the above effects closely follows those in earlier works (e.g. 
YM04; Bromley et al. 2004; Sesana et al. 2004; VR06): we use Monte Carlo merger trees to 
track the hierarchical growth of DM halos, and a simple semi-analytical model to follow the 
growth and dynamics of BHs. We expand over earlier works by adding an explicit calculation 
of the orbits of kicked BHs, and self-consistently include their corresponding time-dependent 
accretion rate. Additionally, we extend the merger tree to redshifts beyond z > 40, and we 
examine a large range of different models. For each set of model parameters, we apply 
our "tree-plus-orbits" algorithm to the entire halo population at z ~ 6 to construct a full 
population of SMBHs at this redshift, and calculate physical quantities of interest: the mass 
function, the SMBH-to-halo mass ratio, the fraction of DM halos hosting SMBHs, and the 
expected detection rate of SMBH mergers by LISA. 

This algorithm assembles SMBHs through simple prescriptions of the aforementioned 
four categories of physical contributions to SMBH formation. We model the seed population 
by assuming that some fraction /seed of DM halos reaching a threshold virial temperature 
Tseed forms a Pop HI remnant BH. In-between mergers, the BHs are assumed to accrete 
gas at a rate of /duty times the Eddington rate mEdd = (1 — e)/e x -^Edd/c^- Here /duty 
should be interpreted as the mean gas accretion rate (averaged over time-scales comparable 
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to the Hubble time) in units of the Eddington rate. Note that this prescription makes no 
distinction between episodic accretion near the Eddington rate during a fraction /duty of the 
time (with no accretion in-between), and constant accretion at all times at a fraction /duty of 
the Eddington rate. We assume that the mergers of BH binaries closely follow the mergers 
of their host halos (but allow for a delay in the latter due to dynamical friction). Finally, we 
simulate the orbits of recoiling BHs under different assumptions about the baryon density 
profile and binary spin orientation. We discuss the relative importance of assembly model 
parameters on the final SMBH mass function and the LISA data stream, and ask whether 
LISA will be able to uniquely determine the underlying assembly model from data. We 
also examine several variants of the above scenario, in which (i) the seed BHs are massive, 
~ 10^ Mq, and formed from the super-Eddington accretion of a collapsed gaseous core; (ii) 
the DM halo is initially devoid of gas when the seed BHs is formed; (iii) seed BHs stop 
forming below some redshift; and (iv) models which maintain the so-called Mbh — o" relation 
between BH mass and (halo) velocity dispersion (Ferrarese & Merritt 2000; Gebhardt et al. 
2000) at all redshifts. 

This paper is organized as follows. In § 2, we detail our methodology by describing 
our assumptions, the assembly algorithm, including the prescriptions of the aforementioned 
physical effects, and the different assembly scenarios we consider. We present and discuss 
our main results in § 3. In § 4, we summarize our most important results, and comment on 
future prospects to understand SMBH assembly at high redshift. To keep our notation as 
simple as possible, throughout this paper the capital M will be used to denote halo masses, 
and m will refer to BH masses. In this paper, we adopt a ACDM cosmology, with the 
parameters inferred by Komatsu et al. (2008) using the five-year data from the Wilkinson 
Microwave Amsotropy Probe {WMAP5 ): ^cdm = 0.233, Qh = 0.046, ^a = 0.721, h = 0.70, 
and (Tg = 0.82. We use Us = I for the scalar index. 



2. Assumptions and Model Description 

Our strategy is as follows: (1) We use Monte Carlo merger trees to construct the hierar- 
chical merger history of DM halos with masses M > IO^Mq at redshift z = 6, i.e. those that 
can host SMBHs of mass m > 10^ Mq; (2) We insert seed BHs of mass ruseed into the tree in 
some fraction /seed of halos that reach a threshold temperature Tgeed! (3) We follow the subse- 
quent BH assembly history by allowing the BHs to grow by gas accretion in-between mergers, 
and by calculating the orbit and accretion history of each recoiling BH in its host halo. We 
assume that BHs add their masses linearly upon merging, and ignore mass losses due to 
gravitational radiation, as these losses never accumulate to significant levels, even through 
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repeated mergers (Menou & Haiman 2004). We prescribe spin distributions of the BHs and 
gas distributions within their host halos. We repeat this procedure for different mass bins 
of the halo mass function, until we have a statistically robust sample to represent the global 
SMBH mass function at redshift z = 6 to an accuracy of a factor of two. We also record 
the BH binary mergers whose masses lie within the mass range ~ (lO'^-lO^) M0/(1 + z), 
corresponding roughly to LISA^s sensitivity range. 



2.1. The Merger Tree 

We construct DM merger trees based on the algorithm by Volonteri et al. (2003), which 
allows only binary mergers. Similar numerical algorithms (e.g. Somerville & Kolatt 1999; 
Cole et al. 2000) give somewhat different results, as have been discussed recently by Zhang 
et al. (2008). We will not reproduce the Volonteri et al. (2003) recipe here in full; instead 
we present a brief review. We take the extended Press-Schechter (EPS) mass function (Press 
& Schechter 1974, Lacey & Cole 1993), 

exp<'4M#^A(f^^, (2) 
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which gives for a parent halo of mass Mq at redshift zq the number of progenitor halos in 
the range M ± dM/2 at a higher redshift z. Here o"m is the amplitude of the linear matter 
fluctuations at redshift z = 0, smoothed by a top hat window function whose scale is such 
that the enclosed mass at the mean density is M (computed using the fitting formula for 
the transfer function provided in Eisenstein & Hu 1999), and 5^ is the redshift-dependent 
critical overdensity for collapse. We take the limit as Az ^ z — zq ^ Q: 

dN _ 1 Mq dal^ , 2 2 N-3/2 d6c , . 

dM~ ^2^M dM ^^^ ^^»^ dz ^'- ^^' 

The two advantages of taking this limit are that (i) by linearizing the expression, the z— 
and M-dependences separate, allowing a tabulation as functions of the parent & progenitor 
halo mass, and (ii) separating A2; allows for a simple algorithm for adaptive timesteps to 
make sure that fragmentations produce binaries at most (no triplets). 

For a parent halo of mass Mp after a small step A2;, the mean number of "minor" 
fragments in the mass range M\o < m < Mp/2 is given by 



M,/2 ^^ 

dM' 



Np= I -jjjdMocAz. (4) 
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We choose Az adaptively such that A^p <^ 1, which ensures that multiple fragmentations are 
unlikely in a given single time step. We place a lower limit of 10"'^ (in redshift units) for the 
timestep to keep computation times manageable. The integral in Equation (jlj) diverges as 
M\o -^ 0, making it computationally prohibitive to compute the merger history of arbitrar- 
ily small halos. To avoid this problem, all progenitors below a fixed mass resolution Mres 
are considered jointly as accreted mass and not tracked individually. Any progenitor with 
M < Mj.es{z) is thus discarded from the tree and its prior history is disregarded. For our cal- 
culations, we choose M^^^iz) to be the mass corresponding to a virial temperature of 1200K, 
corresponding to M^es ~ 4.7 x IO^Mq at z = 40 and 3.4 x IO^Mq aX z = 10. Theoretical 
studies have concluded that Pop III stars can start forming at lower virial temperaturesp 
but numerical considerations have forced us to adopt a somewhat larger threshold. We do 
not impose an explicit upper redshift limit, and we run the tree until our halos at z = 6 
are entirely broken up into progenitors with M < M^^s- As a check on our Monte-Carlo 
merger tree algorithm, in Figure [1], we present the progenitor mass functions of a lO^^M©, 
zq = & parent halo at redshifts of z = 8, 13, 21 and 34, together with the Poisson errors of 
the merger tree output and the predictions from the EPS conditional mass function (eq. 2). 
Our merger tree results are consistent with the EPS conditional mass function up to redshift 
z ~ 40, with agreement within a factor of two for most mass bins and redshift values. In 
particular, the numerical mass function agrees well with the EPS prediction for the low-mass 
progenitors, even at very high redshit, but the higher- mass progenitors are under-predicted 
by a factor of up to two. We note that Cole et al. (2000) used a numerical algorithm similar 
to the one we adopted, to construct a halo merger-tree. As discussed in Zhang et al. (2008), 
that algorithm results in a similar inaccuracy. 



2.2. The Initial Black Hole Population 

The conditions under which the first black holes form are highly uncertain, though nu- 
merical simulations (Abel et al. 2002, O'Shea & Norman 2007) do provide useful indications. 
We parametrize our ignorance in terms of a seeding fraction, such that a fraction /seed of all 
halos reaching the critical virial temperature Tgeed form a seed BH. There are physical mech- 
anisms that make a low seeding fraction plausible: the first stars may form only in rare, 



^Haiman et al. (1996), Tegmark et al. (1997) and Machacek et al. (2001) suggest a threshold virial 
temperature of ^ 400K for collapse. In their recent high-resolution numerical simulations, O'Shea and 
Norman (2007) find star formation in halos of masses (1.5 — 7) x lO^M© between 19 < z < 33, with 
no significant redshift dependence on the mass scatter. These values correspond to virial temperatures of 
260 - 1300K. 



baryon-rich overdense regions with unusually low angular momentum, and seed remnants 
may receive ejecting kicks from collapse asymmetry mechanisms similar to those responsible 
for high-velocity pulsars. Furthermore, radiative and other feedback processes may prohibit 
H2-formation, cooling, and star-formation in the majority of low-mass minihalos at high 
redshift (e.g. Haiman, Rees & Loeb 1997; Mesinger et al. 2006). Since the LISA event rate, 
especially at the earliest epochs, will depend primarily on the abundance of BHs present, it 
is highly sensitive to the seeding function. 

We choose two fiducial seeding models, the first with Tgeed = 1200K (the minimum value 
required for metal-free molecular line cooling and star formation) for a Pop-Ill remnant seed 
BH with mseed = IOOMq. The second model has Tgeed = 1.5 x lO^K and rrisccd = IO^Mq, 
inspired by the "direct collapse" models of more massive BHs from the central gas in halos 
with a deep enough potential to allow atomic cooling (Oh & Haiman 2002; Bromm & Loeb 
2003; Volonteri & Rees 2005; Begelman et al. 2006; Spaans & Silk 2006; Lodato & Natarajan 
2006). If Eddington accretion is the main mode of growth, then we do not expect the choice 
of seed mass for each type of model to qualitatively affect our results, other than the obvious 
linear scaling of the overall BH mass function with the initial seed mass. Only the binary 
mass ratios affect recoil magnitudes, and the subsequent orbital dynamics depends minimally 
on the BH mass. 



2.3. Baryonic and Dark Matter Halo Profiles 

The DM profile for the earliest halos is found to be similar to the NFW (Navarro, Frenk 
and White 1997) profile of lower-redshift, more massive DM halos (Abel et al. 2000; Bromm 
et al. 2002; Yoshida et al. 2003). However, the composition and spatial distribution of the 
baryons, at the time when the seed BH appears in these halos, is poorly understood, and 
is unconstrained by observations. This is unfortunate, since these quantities play a pivotal 
role in determining the orbital dynamics and growth rate of a recoiling BH. 

A steep profile with a cusp will retain BHs more effectively, owing both to a deeper 
gravitational potential well and a larger dynamical drag force at the halo center. The baryon 
distribution will also determine the accretion history of the central BH by determining the 
accretion rate as the BH rests near the halo's potential center, or as it oscillates in a damped 
orbit through the halo following a recoil displacement event. 

In addition, whether the baryons are gaseous or stellar has nontrivial consequences, 
owing to the difference in the dynamical friction force between the two cases. A collisional 
medium provides a greater dynamical friction force than a stellar or DM medium with 
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the same density profile (Ostriker 1999; Escala et al. 2004). Because of tlie difference 
in the drag force, an environment dominated by gas, and not by stars (or dark matter), 
lias several possible consequences on BHs: (1) binaries coalesce more rapidly; (2) a BH 
recoiling in gas has a higher likelihood of being retained in its parent halo; and (3) any 
"vagrant" BH that is displaced from the baryon-rich center of the gravitational potential of 
its host halo takes less time to return there, reducing episodes of suppressed accretion. In 
three-dimensional simulations of star-formation in metal-free minihalos suggest that star- 
formation is inefficient, with either a single star, or at most a few stars, forming at the center 
of the halo (Abel et al. 2000; Bromm et al. 2002; Yoshida et al. 2003, 2008). Since in the 
context of this paper we are concerned with the pre-reionization Universe, we work with the 
assumption that stars are rare before z > 6 and that the baryons in our halos are mostly 
gaseous. 

We model each galaxy as a spherically symmetric mass distribution with two compo- 
nents: a DM halo with an NFW profile, and a superimposed baryonic component. Previous 
studies on this subject (see e.g. Volonteri et al. 2003; Madau & Quataert 2004) have often 
assumed a non-collisional singular isothermal sphere (SIS) profile for the mass distribution. 
This is justified if the gas does not cool significantly below the virial temperature of the 
DM halo, and if it has little angular momentum (so that it is supported thermally, rather 
than by rotation). In most halos whose virial temperature is above lO^K, this assumption 
is less justified, and a disk may form at the core of the DM halo (Oh & Haiman 2002). 
The direct collapse scenario in Begelman et al. (2006) and also VR06 adopt such a "fat 
disk" configuration. However, the central densities of such disks are within the range of our 
adopted spherical profiles. For simplicity, here we only consider three different prescriptions 
for spherical gas distribution. Our fiducial gas profile is a cuspy, p oc r~^'^ power law, where 
we have taken the power-law index of 2.2 as suggested by numerical simulations of the first 
star-forming minihalos (Machacek et al. 2001). This profile is established in halos that 
are able to cool their gas via H2, and describes the gas distribution at the time of the first 
st ar-f ormat ion . 

It is possible, however, that the typical seed BHs are surrounded by a very different gas 
distribution, at the time of their formation. First, the progenitor Pop-Ill stars of the first 
seed BHs are here assumed to form in DM halos of mass ~ 1O^~^M0. The UV radiation from 
the star will photo-heat, and easily blow out most of the gas from low-mass minihalos, even 
before the star collapses to leave behind a seed BH (e.g. Whalen et al. 2004). In this case, 
the remnant BH will find itself in a DM halo devoid of gas, and can only start accreting once 
a merger with another, gas-rich halo has taken place, or until the parent halo has accreted 
enough mass to replenish its gas (e.g. Alvarez et al. 2008). We therefore make the simple 
assumption that no gas is present, until the minihalo containing the newly-formed seed BH 
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merges with another halo, or grows sufficiently - assumed here to be a factor of 10 - in mass. 
However, we will examine the consequences of this assumption below, by performing runs 
without such a blow-out phase. 

Second, as mentioned above, feedback processes may prohibit H2-formation and cooling 
in the majority of the low-mass minihalos (e.g. Haiman, Rees & Loeb 1997; Mesinger et 
al. 2006). The gas in such minihalos remains nearly adiabatic, and can not contract to high 
densities. To allow for this possibility, we will study a variant for the effective gas profile. 
Specifically, we adopt the gas distribution in these halos by the truncated isothermal sphere 
(TIS) profile proposed in Shapiro et al. (1999), which has an r~^ profile at large radii, 
but has a flat core at the center owing to the central gas pressure. The density profile is 
normalized (here, and also in our fiducial gas profile above) such that the total baryon-to-DM 
mass ratio inside the virial radius r2oo equals the cosmological value. Both the DM and the 
baryonic components are assumed to extend out to lOrvir, at which point the density falls 
to the background value. This is consistent with infall-collapse models of Barkana (2004). 



2.4. Mergers of Dark Matter Halos and Black Holes 

We next have to make important assumptions about the treatment of mergers between 
dark matter halos and their resident BHs. 

First, we consider the merger between two DM halos, with the more massive halo referred 
to as the "host" and the less massive as the "satellite" halo. The Press-Schechter formalism 
and our merger tree consider as "merged" two halos that become closely gravitationally 
bound to each other. However, if the mass ratio of a halo pair is large, then in reality the 
smaller halo can end up as a satellite halo, and its central BH will never merge with that of 
the more massive halo. We therefore require in our models that for BHs in such halo pairs 
to coalesce, the halo merger timescale must be shorter than the Hubble time. We take the 
standard parametrization of the merger time: 

r^ Ml r^(^-\ Ml, 

Emerge ^~^ -^ j. <- ^dyn '^ U.iX ^'Yiuhi 

M2 M2 

where M1/M2 > 1 is the ratio of the halo masses, a; ^ 1 is some dimensionless factor that 
encodes the orbit geometry (e.g. circularity) and dynamical friction, and Tdyn and tnub are 
the dynamical and Hubble times, respectively. It has been suggested (Boylan-Kolchin et 
al. 2006) that radial infall along filaments may be preferred in the mergers of elliptical 
galaxies. Boylan-Kolchin et al. (2008) give a fitting formula for the merger time based on 
numerical simulations. Their Equation (5) reduces approximately to Tmcrgcr/Tdyn ~ 0.45 for 
a moderately non-circular orbit with circularity (defined as the ratio of the orbit's specific 
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angular momentum to the angular momentum of the circular orbit with the same energy) 
of 0.5. We take a moderate value of a; ~ 0.5. That is, if M1/M2 < 20, then the BH in the 
smaller halo is considered to be "stuck" out in the orbiting satellite halo and never merges 
with the central BH of the more massive halo. This choice also ensures that the vast majority 
of BH binaries in our simulation do not have extreme mass ratios, as the BH masses co-evolve 
with the host halos. 

We next make assumptions regarding the timescales involving BH dynamics in their 
host halos, as follows: (1) if the two merging halos each contain a central BH, the two holes 
are assumed to form a binary efficiently, i.e. we assume there is no delay, in addition to the 
time taken by the DM halos to complete their merger; (2) the binary is then assumed to 
coalesce in a timely fashion, prior to the interaction with a third hole; and (3) the binary 
coalescence is assumed to take place at the center of the potential of the newly merged 
halo. The ffist assumption has been addressed by Mayer et al. (2007), who report that the 
increased drag force of gas in wet mergers allows the timely formation of supermassive BH 
binaries. The second assumption is valid for binaries in extremely gas-rich environments (see 
Milosavljevic & Merritt 2003 for a review), or in triaxial galaxies (Berczik et al. 2006). As 
for the third assumption: given that the timescales of orbital damping are comparable to the 
intra- merger timescale for BH velocities of > cisiS) unperturbed BHs free-falling during the 
halo merger process are likely to end up near the center of the potential within the merging 
timescales of their hosts. We do not include triple-BH interactions in our analysis. 

The assumption that BH binaries merge efficiently following the mergers of their host 
halos may be unjustified in our models in which initially, the DM halo is devoid of gas, since 
gas is generally believed to be necessary for prompt coalescence. However, this inconsistency 
will not affect our conclusions, for the following reasons. First, we find that a BH merger 
in a gasless environment is a rare event, as it occurs only if both parent halos are low-mass 
halos that had formed seed BHs relatively recently (i.e. neither halos have yet grown in 
mass by a factor of 10). Second, members of such binaries will have equal (or, in actuality, 
similar) masses, since they have not been able to add to their seed mass by accretion. If the 
BHs can merge efficiently without gas, the coalesced product will likely be ejected, owing to 
a shallow halo potential and a relatively large kick of an equal-mass merger. If the binary 
does not coalesce efficiently, it will coalesce once the parent halo merges with a gas-rich, 
BH-free halo, or once the parent halo accretes enough gas to facilitate the merger. Such 
belated mergers will also presumably take place with a mass ratio of close to unity and will 
likely result in ejection, regardless of whether significant gas accretion takes place before 
coalescence. Now, consider the case of a stalled binary encountering a third BH before gas 
enrichment of the halo. If the third BH is much more massive, it will not be ejected by 
gravitational interaction or recoil. There will be a massive BH in the center of the host 
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halo, and whether the two smaller seed BHs were ejected or swallowed by the larger BH is of 
little consequence to our analysis, especially given the rarity of double-seed binaries. Thus, 
inefficient binary coalescence is of concern only when a double-seed binary encounters a third 
BH of comparable mass before the host halo is gas-enriched. Such triple-seed systems are 
likely to be extremely rare indeed, and unlikely to affect the overall mass function at z = Q. 
We anticipate that the main effects of a gas-depleted host halo will be increasing the number 
of similar-mass mergers following the initial epoch of seed formation, and slightly reducing 
the time available for gas accretion. 



2.5. Gravitational Recoil 

2.5.1. Probability Distribution of Kick Velocities 

To calculate the recoil velocities of coalesced BHs, we employ the formulae provided in 
Baker et al. (2008), which are used to fit their numerical results, 

^recoil = ^m + ^i + ^il + 2t^± {Um COS ^ + V\\ siu ^) , (5) 

Vm = Vv^r^4^(l + 5r/), (6) 

2 

VaW = H^{aicos9i- qa2Cos92) , (7) 



3 

-Kjjz [fli sin 6i cos (0i - $1) - qa2 sin 62 cos (02 - $2)] 



where q = m2/mi < 1 is the binary mass ratio, rj = g/(l + qY is the symmetric mass 
ratio, and 9i^2 are the angles between the BH spin vectors 01^2 = 5'i,2/'"^i,2 and the binary 
orbital angular momentum vector. The angles 0i^2 denote the projection of the spin vectors 
onto the orbital plane, measured with respect to a fixed reference angle. As seen from the 
equations themselves, Vm is the kick component that depends only on the symmetric mass 
ratio; f ^ii and v^^ depend on the mass ratio and the projection of the binary spins parallel 
and perpendicular, respectively, to the orbital angular momentum. |f| <l>i(g) = <I>2(1/q') 
are constants for a given value of q that encode the dependence of the kick and orbital 
precession on the initial spin configuration. We use the mean values given in Baker et al. 
(2008) for the fitting parameters: A = 1.35 x 10^ km s-\ B = -1.48, H = 7540 km s~\ 
K = 2.4 X 10^ km s~^, and C, = 215°. We assume spherical symmetry in our host DM halos, 
so we are concerned only with the recoil magnitudes and not with the kick orientations. 



■^Note that this notation differs shghtly from Baker et al. (2008) - we have simphfied their notation to 
be more transparent in our spherically symmetric geometry. 



-13- 



FoUowing Schnittman & Buonanno (2007), for every recoil event, we assign to both 
members of the binary spin magnitudes in the range 0.0 < ai 2 < 0.9, randomly selected 
from a uniform distribution. We consider two scenarios for the spin orientation: a case where 
the orientation is completely random, with < ^12 < vr and < 0i^2 < 27r chosen randomly 
from a uniform distribution; and one where the spins are completely aligned with the orbital 
angular momentum vector o The latter scenario is motivated by Bogdanovic et al. (2007), 
who argued that external torques (such as those provided by a circumbinary accretion flow) 
may help align the binary spins with the orbital angular momentum prior to coalescence, 
making kicks of > 200 km s~^ physically unfavored. While the argument was originally used 
to explain the lack of quasars recoiling along the line of sight at lower redshifts (Bonning 
et al. 2007; although see Komossa et al. 2008 for a candidate recoil detection), the same 
spin-orbit alignment will also impact the early assembly history of SMBHs, by allowing less 
massive halos to retain recoiling BHs. 

Berti & Volonteri (2008) have provided a merger-tree model to follow the spin evolution 
of BHs through gas accretion and binary merger events. In this work, we opt not to track the 
spins of individual BHs due to the uncertainties involved. For instance, if circumbinary disks 
can act to align the spins of each binary member, this would present a scenario significantly 
different from the scenario presented by Berti & Volonteri. As we will show in later sections, 
the spin prescription does not appear to play a primary role in determining the mass function 
oiz = 6 BHs. 

We show the recoil velocity distribution for both orientation scenarios as a function of 
the mass ratio q in Figure [2l The figure shows the mean, l-o", and maximum values for the 
recoil velocity magnitude from 10^ random realizations at a given value of 0.01 < q < 1. 
For q > 0.1, if spins are randomly oriented then kicks for similar-mass mergers are in the 
100 — 1000 km s~^ range, with a handful of kicks above 1000 km s^^ and a maximum 
possible kick of ~ 3000 km s~^; for spins aligned with the orbital angular momentum, kicks 
are typically below 200 km s~^, with the maximum allowed kick no more than 300 km s~^. 
For q ^ 0.1 and for random spin orientations, the maximum kick fmax(Q') is achieved close 
to where v\\ ^ v± is maximized; this occurs when 01^2 = 0.9, cos(0i^2 — "^1,2) = 1, and 
sin6'i ~ — sin6'2 ~ 1. At these spin parameter values, fmax(Q') is well approximated by 



v^ + f y , and is a monotonically increasing function of q. If the spins are aligned with the 
orbital angular momentum vector, then vu = and the maximum kick occurs when ai = 0.0, 



''Both cases have the computational advantage that one does not require the values for $1,2(9)- For a 
totally random spin orientation and a given value of q, choosing 0i_2 randomly is equivalent to choosing 
01,2 — '&i,2 randomly. When the spins are aligned with the orbital angular momentum, $1,2 terms are 
irrelevant because they are always multiplied by sin0i.2 = 0. 
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02 = 0.9. Also, in this aligned case, the spin-independent component Vm is the dominant 
term for q ^ 0.6 and peaks at g ~ 0.345, resulting in the maximum and mean values for the 
recoil speed peaking between these q values. 



2.5.2. Trajectories of Kicked Black Holes 



Given the mass distribution of the host halo and a recoil speed generated from the 
method detailed in the previous subsection, we numerically integrate the equation for the 
radial motion of a BH with mass m, 



dv 

'dt 



GM{r) 



m 

+ aDF -V — , 
m 



(9) 



where r(t) is the radial displacement of the BH from the center of the host halo and v{t) 
is the BH's radial velocity. The first term on the right-hand side is the acceleration due to 
Newtonian gravity with M{r) the total (dark matter + baryon) mass enclosed inside spherical 
radius r; the second is the drag deceleration due to dynamical friction; and the third is the 
deceleration due to mass accretion. A similar calculation of the kicked BH's trajectory has 
been performed by Madau & Quataert (2004) - the main difference from our prescription 
is that they assumed the halo to have a collisionless SIS profile, and adopted parameters 
describing galactic bulges in the local Universe, whereas we use the hybrid DM+gas profile 
described above, and adopt parameters relevant to low-mass halos at high redshifts. 

For a non-collisional medium (in our case, for dark matter), the dynamical friction is 
described by the standard Chandrasekhar formula (see e.g. Binney & Tremaine 1988), 



at^{r,v) 



-A'KG'^mpir) - In A 

V 



erf(X)-^Xexp(-X^ 

7r 



(10) 



where InA is the Coulomb logarithm and X = ^/(v^o'dm), with ctdm the velocity dispersion 
of the DM halo. In a collisional medium, the density wave in the wake of an object traveling 
at near or above the sound speed is enhanced via resonance, an effect that has no counterpart 
in collisionless media. This results in an enhancement of the dynamical friction force, for 
which Ostriker (1999) has derived an analytical formula. However, the Ostriker prescription 
is known to overpredict the drag force at slightly supersonic velocities when compared with 
numerical simulations. While Escala et al. (2004) provides a fitting formula that better fits 
the numerical results at low speeds, their formula suffers from the opposite problem, and over- 
predicts the drag for highly supersonic motion. We therefore adopt a hybrid prescription, 
adopting the Escala et al. (2004) formula for motion with Ai < Aieq and the Ostriker formula 
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for Ai > Aieq, where Aieq is the Mach number (= |f |/cs) where the two prescriptions predict 
equal drag. The entire prescription is described by 



fiM) 



gas 

"df' 

0.5 In A 

1.51nA 



r,v) = -4:7[G'^mp{r)- x f{M), 



erf(^ 



iln(#hi)-lnA 



.M 



^Mexp{-My2) 
^Mexp{-My2) 



where : (11) 

if < A^ < 0.8; 

a 0.8 < M < Meg] (12) 

if M > Meg- 



The Coulomb logarithm In A is not a precisely known parameter, but is generally agreed to 
be > 1 for both the stellar and the gaseous cases. We adopt the value InA = 3.1 used in 
Escala et al. (2004), which yields Meg ~ 1.5. 

The drag force depends on the local gas sound speed. Instead of attempting to compute 
a temperature profile explicitly, we make the approximation that the gas sound speed is 
constant and given by the isothermal sound speed of the halo virial temperature. We believe 
this to be justified from numerical simulation results that show the gas temperature to vary 
by at most a factor of ~ 3 within the virial radius despite a steeper-than-isothermal (oc r~^'^) 
density profile (see, e.g. Machacek et al. 2001). While local variations in the sound speed 
may have significant effects when f ~ c^, the recoil events of interest here are for the most 
part highly supersonic. Recoil events with w ~ Cg will result in quick damping of the BH 
orbit and for the purposes of this paper will in all likelihood be indistinguishable from the 
zero-recoil calculation in terms of their accretion history. 

The virial temperature is given by (e.g. Barkana & Loeb 2001) 



T,ir ^ 370(1 + z)(^ 



0.6 



M 



llO^M^ 



2/3 



0.14 



1/3 



K, 
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and the isothermal sound speed is 



1.8(1 + ^)^/2 



/ M 
VlO^Mo 



1/3 



0.14 



1/6 



km s 



-1 



(14) 



We also employ a simplified prescription for the velocity dispersion of non-collisional matter, 
given by the SIS expression evaluated at the virial radius: ctsis = \/GM/2r2m- The simplified 
expression agrees with the exact velocity dispersion for the NFW profile within ~ 20% inside 
the virial radius. 

Because matter that is bound to the BH does not contribute to dynamical friction, 
we follow Madau & Quataert (2004) and truncate the density profiles at the BH radius of 
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influence, tbh ~ G'm/agjg. The density is furthermore assumed to be constant inside this 
radius. Although the BH will drag with it the surrounding gravitationally bound matter, 
effectively increasing its mass, the additional mass is small owing to the large initial recoil 
velocity (e.g. Lippai, Frei & Haiman 2008), and we have ignored this mass-enhancement 
here. 

Figure [3] illustrates the effect of the halo matter distribution on the BH orbit. For each 
of the three orbits shown in the Figure, we adopt M = IO^Mq, m = 10^ Mq, z = 20 and 
■^kick = 100 km s~^. The black curve shows, for reference, the radial orbit for a pure NFW 
profile. The red curve corresponds to the case which includes an NFW DM component and 
a power-law gas profile with p ex r~^'^. The blue curve is for a halo with an NFW DM 
component and a TIS gas profile. 

In our calculations, we are mainly interested in whether the kicked BH is ejected and 
lost (i.e., can not contribute to the final SMBH mass at z = 6), or is retained (i.e., eventually 
returns to the nucleus, and can be incorporated into the z = 6 SMBH). We place the following 
retention condition for recoiling BHs: the BH must return to within 1/lOth of the virial radius 
of the newly merged host halo within 1/lOth of the Hubble time. The fate of a (SM)BH 
placed in an orbit extending to the outskirts of its host halo is uncertain: even if it is not lost 
through tidal interactions with an incoming merging halo, it is not likely to form a binary 
that hardens efficiently. We therefore impose the above conservative cutoff, in order to avoid 
tracking these vagrant BHs. Our retention threshold velocity, fret? above which recoiling 
BHs do not return within the prescribed time limit and are considered lost, is a numerically 
calculable function of M, m, z and halo composition. In order to minimize computation time, 
we tabulate v^et in the range 5 < z < 40, IO^Mq < M < IO^^Mq and 10"^ < m/M < 1 
and approximate the result with a fitting formula that accurately reproduces the exact 
numerical results within 5% in the tabulated range. In principle, we should compute the 
retention velocity as a function of the time to the next merger experienced by the halo. 
However, we find the time dependence to be weak. The distribution of the time intervals 
between halo mergers in a given merger tree has a sharp peak at ~ lO'^tnubj with far fewer 
mergers occurring at ~ 10~^tHub and ~ tnub- At these tails of the distribution, v^-et varies by 
^ 10% from the value for 10~^tHub- We find that v^^t ~ 5 — 8 x asis- This is comparable to 
the escape velocity for a non-dissipative pure SIS profile that is truncated at the BH radius 
of influence, v^^c ~ Sasis if m = lO'^M (YM04). 

The weak dependence on fret on the return time limit is a counterintuitive result, but it 
can be understood as follows. There is a minimum kick speed that is required to displace the 
hole beyond 0.1r2oo, which represents limi^o'^^ret; and there is a maximum kick, limf^oo'^ret, 
beyond which the BH remains completely unbound from the halo, even in the presence of 
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drag. frct(^), then is a function of time that is always between these two extreme values. 
However, owing to the high central density of our gas-dominated halos, the difference be- 
tween these two limits is small, ~ 10%. Since this difference is smaller than other model 
uncertainties, such as those stemming from discrepancies from the actual density profile of 
high-redshift DM halos (e.g. clumping or triaxiality) and the merger tree prescription, we 
simply use the retention velocity computed for the approximate median time limit, O.ltHub- 

Finally, for simplicity, we treat the halo as a static mass distribution during each recoil 
event. That is, we ignore the cosmological evolution of the DM potential, and we assume 
that the recoiling BH does not affect the medium by clumping or heating it. Note, however 
that the latter feedback may play a nontrivial role in real systems, since the kinetic energy 
of a recoiling hole can be comparable to the gravitational binding energy of the entire host 
halo and can be expected to cause significant disruption of the surrounding matter. 



2.6. The Black Hole Accretion Rate 

We turn now to our prescription for the accretion rate of the BHs in our model. Of 
particular interest is the possibility that the gravitational recoil effect will significantly limit 
the ability of kicked BHs to accrete gas, by displacing them into low-density regions for 
prolonged periods, and/or by limiting through high relative velocities the amount of gas 
that can be gravitationally captured. One can imagine a scenario in which a SMBH whose 
progenitors have survived numerous kicks but have spent long episodes in underdense regions 
may have a final mass much less than that predicted by simple Eddington growth. We 
therefore follow the accretion rate self-consistently, as the recoiling holes proceed along their 
radial orbits. Specifically, in our models a BH embedded in gas is assumed to undergo 
standard Bondi-Hoyle-Littleton (BHL) accretion. 



A-KG'^Pb{r, 2 



' ' ' (c2 + t;2)3/2 

The accretion rate is capped at the Eddington rate. 



m\ (15) 



m = , (16) 

e tEdd 

where tsdd = 44Myr and e is the radiative efficiency, for which we assume e = 0.10. 

If the gas density is too low, or if the sound speed or its velocity with respect to the gas 
disk too high, a BH may not be able to accrete at the Eddington rate even if it is close to the 
center of a halo. Because the BHL rate is proportional to m^, an underfed BH with initial 
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mass mo will eventually reach the Eddington accretion rate (oc m) at a threshold mass 

where 4 km s~^ is the isothermal sound speed for an ionized hydrogen gas at 1200K. The 
typical central density of a 1200K halo with a TIS gas profile is ~ 5 x 10~^(1 + z)^/'^M0pc"^. 
Sub-Eddington accretion rates are not an issue for BHs in a power-law gas profile, as the 
steep profile provides a sufficient central density for immediate Eddington accretion. 

The time it takes for a BIf with tuq < niEdd to reach this threshold mass, assuming that 
the local density remains constant, is 



''crit 



(^ - T^'- - - (i^)"" ii^r (ivi^)"'«- ^- 



tcrit ~ a few 100 Myr for IOOMq halos embedded in 1200K TIS halos at z > 20, and 
^crit ^ Gyr for z ^ 14. This means that for TIS gas profiles, seed holes will spend a 
significant fraction or all of the available time prior to z ~ 6 accreting below the Eddington 
rate. 

The difference between BHL and Eddington accretion rates as they relate to BH growth 
is also discussed in Volonteri & Rees (2006). However, in that paper the context is for 
the direct formation of m > IO^Mq intermediate-mass BHs through super-Eddington BHL 
accretion. We here adopt the opposite extreme assumption, i.e. that the BH radiates 
efficiently at all times, and its accretion rate obeys the Eddington limit. The BHL rate can 
then initially be sub-Eddington in TIS halos, owing to the low BH mass and low gas density 
(the baryon density required to fuel BHL accretion at the Eddington rate was also discussed 
by Turner 1991). 

To illustrate the impact of extended sub-Eddington growth phases, we perform a simple 
analytical calculation. Suppose that a BH is formed with mass mseed at redshift z in a halo 
with virial temperature 1200K, and that the local gas density is held constant at the value 
when the BH was formed. Note that for this exercise, we assume that gas density is constant 
even as the halo around the BH is growing by merging with other halos - in other words, 
we assume these mergers deliver gas to the nucleus containing the original BH, roughly 
maintaining a constant density at the Bondi radius around the BH. Figure H] shows the 
maximum possible mass that can be attained by such a BH growing in isolation through 
gas accretion alone, assuming that the accretion rate is determined solely by Equations (fT5l) 
and (ITB]1 and that accretion is not supply-limited. If the host halo has a steep power-law 
cusp, the accretion rate is Eddington throughout regardless of when the seed BH is formed. 
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However, if the central fuel density is low, then it is possible for the local BHL accretion 
rate to be significantly sub-Eddington initially. In such a scenario, the earliest forming seeds 
are the only ones able to reach the Eddington rate; the late-forming seeds are unable to 
reach the Eddington rate before z = 6. In this scenario the late-forming seeds, which are 
easily identifiable by the drastically shallower growth slope in the figure, cannot grow rapidly 
enough to contribute to the SMBH population. Note that assuming a constant gas density 
in this calculation gives an optimistic accretion rate for the TIS case, as in those profiles the 
central gas density generally decreases with Hubble expansion and significantly reduce the 
BHL rate. 

It is computationally expensive to numerically integrate the individual orbits and accre- 
tion histories of every recoiling BH in our simulations. We therefore tabulate the accretion 
growth in Eddington units during the first O.ltnub of the orbit, and describe the results in 
a fitting formula in the same manner as we have done for the retention velocity. Given a 
specific prescription for the baryon distribution, we tabulate across four relevant variables 
in the following ranges: IO^Mq < M < IO^^Mq, 10"^ < m/M < 1, 5 < ;z < 40, and 

< t^kick/^ret < 1- 

In the absence of a kick, and if the accretion rate were always at the Eddington limit, 
the SMBH mass in a given halo at 2; ~ 6 is easily approximated by 



AT I £ ^ ''seed, 6 \ 

V "^ 1 - e tEdd / 
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where mgeed is the seed BH mass, A'seed is the number of seeds in the merger history of 
the halo, tseed,6 is the available time between the typical seed formation time and z ^ 6 
and /duty is the time-averaged duty cycle for accretion. Equation (TT^ represents the ideal, 
maximally efficient scenario for SMBH assembly, and we can use it to effectively measure 
the cumulative impact of underfed accretion, recoil-induced ejection, and other factors that 
limit the assembly efficiency. 

The measurements of clustering of quasars in the SDSS suggest that the duty cycle of 
active (luminous) accretion increases steeply with redshift at 3 ^ 2; ^ 6, with the most active 
quasar BHs at z ^ Q showing 0.6 < /duty < 0.9 (Shen et al. 2007; Shankar et al. 2008). We 
therefore adopt duty cycles of > 0.6. Although it is likely that SMBHs regulate their own 
growth through feedback mechanisms, we do not address such scenarios a priori. To keep 
our models as simple as possible, we will only impose the loosest upper limit on the SMBH 
accretion rate: they must not accrete more mass than there the total mass of baryons in 
the host halo. However, we will discuss an alternative ad-hoc model below, which is able 
to reproduce the well-known relation between SMBHs and the velocity dispersions of the 
bulges of host galaxies (the m-a relation). 
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2.7. Putting Together the z = 6 SMBH Mass Function 

Explicitly constructing the SMBH mass function aX z = 6 over a wide mass range is 
computationally intractable. The host halo mass inferred from the observed quasar space 
density from the z ~ 6 quasars is several xIO^^Mq (e.g. Fan 2006). For every halo with 
this mass, there are ~ 10^ halos with 10® — 10^ Mq. Blindly calculating the SMBH mass for 
every halo with M > IO^Mq using our trees-plus-orbits algorithm would be prohibitively 
expensive computationally. 

We therefore carry out a piecewise calculation of the SMBH mass function that is 
computationally tractable. The procedure is as follows: (1) We group the halo population 
into logarithmic mass bins of size x < log^^Q M < x+Ax; (2) For each bin, we select > 10^—10'^ 
individual Monte- Carlo- generated halos with masses randomly generated from the Press- 
Schechter distribution at z = 6; (3) We simulate the BH population for each such halo using 
our trees-plus-orbits algorithm, and assume the resulting sample is representative of all ;i; = 6 
halos in the mass bin; (4) For each bin we multiply the sample by the appropriate weight to 
construct the entire Press-Schechter halo mass function, J dN/d In Md In M; and finally 
(5) Sum the contributions from each bin. The result is a Press-Schechter distribution of 
halos with M > IO^Mq aX z = 6, with a statistical representation of the corresponding 
SMBH mass function. The bins used and their relevant properties, including the number of 
Monte-Carlo halos that were cloned to populate the full mass function, are listed in Table 1. 
This numerical shortcut is not used for the most massive halos. 40 halos are expected above 
IQ^'^-^^Mq, and these are simulated individually. 

This method allows a fast calculation of the BH mass function, with the caveat that 
there must be enough BHs in each bin to keep statistical uncertainties to a minimum. Unfor- 
tunately, this is not always preventable for models with extremely low /sccdj and the reader 
will notice statistical noise in the results of such models. In some cases, we increase the halo 
sample by a factor of 10 in an attempt to reduce the errors. 

In all, our simulations represent the Press-Schechter population of DM halos in a co- 
moving volume of ~ 280Gpc^, roughly equal to the co moving volume that was probed by 
the SDSS between z ~ 5.7 — 6.4. We have simulated a total of ~ 1.17 x 10^ DM halos (see 
the column Ngim in Table 1), and through the procedure described above extrapolated the 
results to represent the ?a 3 x 10^^ Press-Schechter halos (with M > IO^Mq) expected to be 
present in the 280Gpc^ volume. 
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3. Results and Discussion 

3.1. Building the > IO^Mq SMBHs 

As stated in the Introduction, our primary goals are (i) to understand the possible 
ways in which the > lO^M© SMBHs may have been assembled at redshift z > 6, and (ii) 
whether the LISA event rates are sufficiently different in competing models so that one 
can disentangle the physical assembly scenario from the LISA data stream. As a ffist step 
toward these goals, we would like to survey all feasible combinations of the physical assembly 
parameters, understand the impact of each model parameter, and look for corresponding 
give-away features in the predicted LISA stream. 

Although a broad simulation survey is beyond the scope of this paper, we have un- 
dertaken a cursory tour of the basic parameter space. We begin by considering two basic 
seed models: IOOMq seeds forming as remnants of Pop III stars in minihalos when they 
first reach the virial temperature T^u > Tseed = 1200K and IO^Mq seeds forming as a result 
of direct collapse of gas in halos when they first reach Tvir > 1.5 x lO'^K. For each case, 
we consider halos with an NFW DM component and a gaseous component with either a 
cuspy Pgas oc r~^'^ power-law or a corey TIS profile. The TIS models consistently failed to 
produce SMBHs hj z = 6, with typical maximum BH masses of ^ 10^ Mq. In those models 
the central gas densities are too low to allow for prolonged episodes of accretion near the 
Eddington rate, as noted in Section 2.5 above. 

For each type of seed, we vary three sets of parameters: (i) the seeding fraction 10^^ < 
/seed < 1, i-e. the probability that a halo reaching the critical temperature will form a seed 
BH; (ii) the time-averaged accretion rate in Eddington units, characterized by the duty cycle 
/duty, for which we use /duty = 0.6, 0.8 and 1.0 (note that /duty and the radiative efficiency e 
are degenerate in our prescription; see below); and (iii) whether at the time of their merger, 
the BH spins are randomly oriented or aligned with the orbital angular momentum of the 
binary. The spin magnitudes are chosen from a uniform random distribution < ai^2 < 0.9 
(Schnittman & Buonanno 2007). We do not track the evolution of BH spins in our models. 
While more rapidly spinning BHs are capable of accreting at higher efficiency, we neglect this 
effect and assume a global efficiency coefficient in our models. Of the four main ingredients 
of the SMBH assembly introduced in § 1 above, we here therefore vary three: the seed rarity, 
the accretion rate, and the recoil dynamics. For the fourth, the merger rate, we have simply 
assumed that BHs merge when their parent halos do, as extreme-mass BH binaries are rare 
in our simulation, given the threshold we have imposed on the mass ratio for halo mergers. 
We will call the above our fiducial set of parameters. 



For each realization, we compute the mass function of BHs at z = 6 with m > lO^Mg 
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the rate of events in the LISA observational mass range, ~ 10^~^/(1 + z)Mq; the fraction 
of DM halos hosting a central massive BH; and m/M, the ratio between the mass of the 
SMBH and its host DM halo, which serves as a proxy for the m-a relation. 

We begin with Figure [5], showing the mass function for the mseed = IOOMq Pop III 
seed model. This and the companion figures are organized with different values for /seed in 
different columns, the two spin prescriptions in different rows, and the different duty cycles 
in different line styles (and different colors, in the online version). This will be the default 
organization of our 8-panel figures unless noted otherwise. 

The numbers in the upper-right-hand corner of each panel represents the total SMBH 
mass density, \ogiQ[p,/{MQ Mpc~^)], for all BHs with m > IO^Mq. Because of statistical 
fluctuations, for multiple model realizations with identical parameters this value can vary by 
^ 10% for /seed ^ 10~^, and as much by a factor of a few for /seed ^ 10~^. For each of these 
simple models, this density is exceedingly high compared to the SMBH density of the local 
universe. For the present discussion, we will overlook this point as we address the effects of 
the various model factors; we will introduce the additional constraint from p, in subsequent 
sections. 

The most stringent observational requirement for the high-mass end of the z ^ 6 SMBH 
mass function comes from the SDSS observation of the z ^ 6.4 quasar J1054+1024, which 
has an inferred mass of ~ 4 x IO^Mq. Since this object was detected in a region ~ 10% of 
the sky, we estimate that > 10 similar objects may exist at z ~ 6. We have adopted this 
as a rough indication of the lower limit of the SMBH mass function at redshift z = 6, and 
represent this limit by the upper right quadrangle in each of the panels, delineated by the 
red dashed lines. Note that these lower limits are conservative, since they do not require 
the presence of any additional SMBHs with comparable mass that are "off". Since high 
values for the duty cycle - near unity - are suggested by quasar clustering measurements 
(as discussed above), and are also required for growing the most massive SMBHs (as we find 
below), this correction is only of order a factor of ~two. As seen in the figure, the high-mass 
end of the SMBH mass function is so steep that increasing the lower limit on the required 
SMBH space density by even ~ 2 orders of magnitude would make little difference to our 
conclusions. 

The first conclusion one can draw from this figure is that several combinations of pa- 
rameters can produce SMBHs massive enough to power the quasar J1054+1024. We will 
return to this and other observational constraints in the next subsection, but let us for now 
focus on understanding the effects of our various parameters and their combinations. In gen- 
eral, the effect of varying each of the parameters is relatively straightforward to understand. 
Increasing the accretion rate, increasing the seed occupation fraction, and aligning the spins 
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all tend to result in more massive BHs. Note that the accretion rate in Eddington units, 
/duty, is degenerate with the accretion efficiency, as m oc /duty >< (1 ~ e)/e- For example, a 
model with e = 0.15 and /duty = 1-0 is equivalent to e = 0.10 and /duty = 0.63. We have 
used e = 0.10 throughout the results presented here. With this value for the efficiency, it 
is possible to build the most massive SDSS quasar SMBHs by 2; ~ 6, starting with lOOM© 
seeds. If the efficiency is ~ 0.15, however, it is only possible to build the SMBHs in question 
with the most optimistic assumptions: /duty ~ 1, /seed ^0.1, and spin alignment would all 
be required. 

We note two non-trivial observations to be made from Figure [51 First, if the seed 
fraction is low, the spin orientation has a minimal effect on the BH mass function. This 
is because if seeds are extremely rare, they are likely to grow in isolation for much of their 
existence along with their host halos. The ffist mergers are not likely to occur until the 
gravitational potentials of their host halos are deep enough to retain them from any recoil 
event. Second, the increase in the SMBH abundance from increasing the seeding fraction 
has a tendency to plateau. This is the reverse situation compared to the low /seed limit: if 
seeds are very common, they are likely to experience multiple BH mergers very early in the 
merger tree, when their masses are still comparable, and ejections become very common. 
As /seed iucreascs, assembly becomes increasingly inefficient at early times. In models where 
the seeding halo temperatures are lower, the efficiency saturates at lower values of /seed- 
Furthermore, in models with the shallower TIS gas profiles, we find that increasing the 
occupation fraction beyond a certain "sweet spot" value rapidly decreases the final SMBH 
masses. The reasons for this are that (i) the seed BHs in these models barely grow, making 
near-equal mass (g ~ 1) BH mergers, and therefore large kicks, much more common and 
(ii) the gas drag is reduced, making it easier to kick holes out of these halos. We conclude 
that arbitrarily increasing the number of seed holes contributing to the assembly process 
is not necessarily an efficient solution to the SMBH assembly problem. In fact, excessive 
seeding can lead to a different conflict with observations - overproducing the mass density 
in lower-mass SMBHs -, to which we will return in the next subsection. 

Let us turn now to Figure [6], which shows the BH occupation fraction in M > 10^ Mq 
DM halos at 2 = 6, and the BH-to-halo mass ratio (which here serves as a proxy for the m-a 
relation; see, e.g. Ferrarese 2002) for each of the models. Here, we have arbitrarily defined our 
occupation fraction to mean the fraction of DM halos that host a BH with a minimum mass 
of IO^Mq, as we are interested in the population of nuclear BHs and not stellar remnants. 
Menou et al. (2001) have shown that, in the low-redshift Universe, the fraction of DM halos 
hosting a SMBH will approach unity, regardless of the initial occupation fraction at earlier 
times. We find that at z ~ 6, the occupation fraction in halos with M ^ IO^^Mq can still 
be significantly below unity, depending primarily on /seed- In principle, a future survey that 
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can determine the quasar luminosity function (LF) at 2; = 6 to several magnitudes deeper 
than the SDSS could look for this drop in the occupation fraction, since it will produce 
a flattening of the LF at magnitudes below some threshold. The spin prescription has a 
noticeable effect on the z = Q occupation fraction for /seed ^ 10~^. On the other hand, the 
duty cycle essentially only affects the mass of the BHs, and not their presence or absence, and 
so has a minimal effect on the occupation fraction, within the parameter range shown. Note 
that we do not expect to reproduce the m-a relation in our simulations, as we employ simple 
and ^-independent accretion and seeding prescriptions, and our models have no feedback to 
enforce the relation. Instead, the m/M relation should be taken as a sanity check that we 
are producing a physically viable BH population. Note that in some of the models shown in 
Figure 5, the BHs grow much larger than the m/M relation observed in nearby galaxies. In 
particular, as the figure reveals, SMBHs in the lower-mass (~ 10* Mq) halos in the models 
with /duty = 1 tend to consume most of the gas in their parent halos, clearly an improbable 
result. 

For an explicit comparison, alongside the m/M relation produced by our models, we 
have plotted the the value expected based on measurements of the m,-a relation in local 
galaxies. Specifically, we show the expression from Wyithe & Loeb (2003), 

We adopt their parameter choices of eg = lO"^"^ and 7 = 5. This expression satisfies the SDSS 
constraints at the high-mass end of the mass function. It also agrees well with the relation 
suggested by Ferrarese (2002) for SMBHs in the local universe, m ~ 1O'^(M/1O^^M0)^-^^. As 
the figure shows, our predicted m/M relation tends to have the opposite slope than the one 
inferred from the observed m-a relation: in our results m/M decreases with mass or stays 
roughly constant as M increases, but this is the opposite of the empirical trend. This is due 
mainly to the difference in the growth rates of halos and holes: our simple prescription for 
steady, exponential accretion for the BHs can significantly exceed the growth rate of DM halo 
masses due to the accretion of unresolved low-mass halos in the EPS merger tree. As a result, 
in some cases, the host halo mass predicted for the z = Q quasars is as low as IO^^Mq, which 
is an order of magnitude lower than would be predicted from the extrapolation of the locally 
measured m/M relation, or from the inferred space density of the host halos (e.g. Haiman & 
Loeb 2001). However, any extrapolation of the m — a 01 m/M relation to high redshift, and 
the masses of halos that host the brightest z > Q quasars, at present, have large uncertainties, 
and do not robustly preclude such low values. As will be discussed below, the overly rapid 
growth of SMBHs in this suite of models motivates modifications to the modeling, including 
a model in which the extrapolated m/M relation holds by assumption. 
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Figure [7] shows estimates for the LISA event rate, calculated for all binary mergers 
satisfying IO^Mq < (mi + m2)(l + z) < 10'^ Mq as (see, e.g. Menou et al. 2001), 

d'^N , ,2 , , AiV 

^j^ = 4vrc4o.W^^^, (21) 

where AA^ is the number of SMBH merger events in the tree in a time step Az and a 
simulated comoving volume AV, and dcom is the comoving distance. Although there is a 
mild dependence on the duty cycle / accretion rate and the kick prescription, it is evident 
that for our assembly models /seed has the greatest effect in setting the rate of SMBH binary 
mergers detectable by LISA. Because the initial merger rates scale as fseed^ ^^^ measured 
event rate is extremely sensitive to the BH number population. Since the merger activity 
peaks near 2; ^ 10, LISA should be able to measure the masses of most of these SMBH 
binaries to high prevision (Hughes 2002, Arun et al. 2008). Although the raw detection rate 
- without any information on BH spin or mass ratio - alone will be richly informative, the 
ability to determine the source masses with high fidelity should be very useful in further 
constraining the growth rate of the first SMBHs, and breaking degeneracies between the 
allowed parameter-combinations. The observed distribution of binary masses as a function 
of redshift will provide direct snapshots of the mass function and shed light on its evolution, 
independently from quasar luminosity surveys (e.g. Yu & Tremaine 2002; Willott et al. 
2006). 

Another family of assembly models that has been frequently discussed in the literature 
is the so-called "direct collapse" model, wherein BHs with m ~ 1O^~^M0 are formed rapidly 
from gas cooling via atomic H in halos with virial temperature T > lO^K (Oh & Haiman 
2002; Bromm & Loeb 2003; Volonteri & Rees 2005; Begelman et al. 2006; Spaans & Silk 
2006; Lodato & Natarajan 2006). We simulate such a family of models, for the same set of 
the parameters /seed and /duty and the same spin alignments. We choose mseed = lO^M© and 
Tseed = 1-5 X lO^K. We show the mass functions, the hole-halo relations and the LISA rates 
(the same information as in Figures O [6] and [7]) in Figures [H [9] and [TOl Although the main 
differences are all fairly intuitive, it is instructive to address how the direct-collapse models 
differ from the Pop-HI seed models. 

First, from Figure [8] we see that it is much easier to construct more massive SMBHs 
owing to the larger seed masses. In fact, the models with /seed ^ 0.1 become unphysical, since 
the SMBHs in these models would exceed the baryon mass {Qb/^m)M of their parent halos. 
The second point is that the mass function is very differently distributed between the Pop- 
III and direct-collapse scenarios. The reader will note that for each set of parameters, the 
overall SMBH density does not differ significantly between the corresponding seed models. 
However, this is deceiving as the mass function is clearly different, with the direct-collapse 
seeds resulting in a more "top heavy" SMBH population. For the range of parameters 
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surveyed, the Pop-Ill model has > 90% of the SMBH mass in the m < IO'^^Mq range if 
D = 0.6, compared to ^ 30% for the D = 0.6 direct-collapse seed models. For D = 1.0 
and random spin alignment, ~ 2 — 5% of the SMBH density resides in the most massive 
m > IO^Mq BHs if the seeds are Pop III; for the same parameters, ~ 60 — 70% of the mass 
is in the billion-plus solar mass BHs in the corresponding direct-collapse models. Note that 
the total BH mass density remains extremely high; again, we will address this point shortly. 

Third, there is a slightly weaker dependence on the spin orientation of the BH binaries. 
This is most apparent by comparing the z = 6 occupation fractions in Figures El and [HI 
and is due to the deeper potentials of the host halos in the direct collapse case. Fourth, 
even though the m/M relation continues to have a slope opposite to the locally observed 
trend, there appears to be a break in the relation at M ~ 10^ Mq, accompanied by a 
drop in the occupation fraction. This is due to the simple fact that halos below this mass 
threshold cannot have many T > 1.5 x lO^K progenitors and the model similarly prohibits 
intermediate-mass BHs. This cutoff contributes to the top-heaviness of the mass function for 
these models. Fifth, LISA rates are lower by one to two orders of magnitude than in the Pop- 
III seed models, because there are fewer 1.5 x lO^K halos than 1200K halos (another factor 
is that the seed BHs are already born with a mass near the middle of LISA's logarithmic 
mass range, so they spend only ~ half the time in the LISA band, compared to the Pop-Ill 
seeds). It is worth noting, in particular, that it is possible to build the SDSS quasar BHs 
in ways that produce no detectable GW events for observation with LISA beyond 2 ~ 6 (in 
contrast, in the successful pop-Ill models, a minimum of a few events are predicted). In 
such scenarios, SMBHs are extremely rare until z ~ 6, at which point the SMBH occupation 
fraction will evolve toward unity fairly rapidly, as described by Menou et al. (2001). 

Finally, in Figure [11] we present the mass function and the LISA detection rate in five 
variants of a fiducial /duty = 0.65, /seed = 1, aligned spin model with pop-Ill BH remnant 
seeds. We choose these values because they just barely are able to match the abundance of 
the most massive SDSS quasar BHs, and because they produce the most BH mergers, and 
thus the effects of varying the other recoil-related parameters will be the most discernible. 
We show the fiducial model in bold. In the modified models, we (i) allow the gas density 
profile to be shallower, with a r~^ power law - this is to allow for the possibility that the gas 
surrounding the BH has not cooled and condensed to high density (dark blue curves); (ii) 
require the BH spins to be near-maximal at ai,2 = 0.9, instead of choosing them randomly 
from the range 0.0 to 0.9 - this is to allow for the possibility that all SMBHs at high z are 
rapidly spinning, e.g. due to coherent accretion (green curves); (iii) allow halos of all mass 
ratios to merge, where previously we had considered a halo with mass less than l/20th that 
of its merger companion to become a satellite (yellow curves); (iv) assume that the Pop III 
seed progenitors are not able to blow out the gas in the host minihalo (red curves); and (v) 
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ignore the effects of accretion suppression due to episodes of recoil-induced wandering, and 
instead assume tliat all BHs accrete at /duty unless ejected (light blue curves). 

We also ran models with a TIS profile for the gas component. We found that these 
models always failed to produce any SMBHs above IO^Mq by redshift z = 6 due to the 
initial phase of sub-Eddington accretion of seed BHs, and their results are not shown in 
Figure [TTJ This implies that SMBHs must be continuously surrounded by dense cores of gas 
that was able to cool at the centers of DM halos - feeding holes with the low-density gas 
in DM halos whose gas was unable to cool does not allow for high enough accretion rates 
(Turner 1991 emphasized the same issue for the growth of 2; ~ 4 quasar BHs). 

The results shown in Figure [TT] give insight to the importance of the assumptions that 
went into our models. In particular, of all the parameters that we have varied for fixed /seed 
and /duty; the most significant for the SMBH mass function, by far, are the spin orientation 
and the limit on the halo mass ratio for timely mergers. Both of these have a similar effect 
of increasing the number of BHs, especially at the massive end. The former result - that 
maximizing the spin increases the SMBH mass function - may seem surprising at first, since 
generally large spins imply larger kicks. However, this is not always the case, as can be 
understood from equations (5-8). Under the assumption that both spins are aligned with 
the orbital angular momentum, i;|| = 0, and v± oc (ai — 502) is maximized for unequal spins 
(i.e. ai = 1 and 02 = for a typical g ~ 1); setting ai = 02 = 0.9 therefore eliminates the 
largest kicks, and allows more BHs to be retained. Comparatively smaller differences are 
visible in the mass function for the other model variants. Figure [TT] also shows that adding 
in the unequal-mass halo mergers and increasing the spins affect the LISA event rates 
differently: the former adds new events mostly at z < 10, where unequal-mass mergers are 
more common, whereas increasing the spin mostly adds new events at z > 10. Interestingly, 
ignoring the blow-out has little impact on the SMBH mass function at z = 6, but it does 
shift the LISA events to higher redshifts, especially at 2; > 10. Figure [TT] suggests that the 
LISA event rate can be useful in disentangling these three effects. 

Perhaps the most surprising inconsequential variation is ignoring episodes of reduced 
accretion due to the BH wandering in low-density outskirts following recoil. The reason for 
this is simply that lengthy oscillating orbits are relatively rare if the central gas density is 
high; the BHs either return quickly, or are ejected (by assumption). Recall that our definition 
for a "retained" BH is that the recoiled hole must return to within 1/lOth of the halo's virial 
radius within 1/lOth of the Hubble time. For low kick speeds, the BH does not get very far, 
because the gas provides both a steep gravitational potential barrier and a strong dissipative 
sink. The orbit will thus be rapidly damped, with only very brief periods of underfeeding. If 
a BH is kicked hard enough to reach the outskirts of the halo, there is very little dissipation 
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there, to tug it toward the center or to further damp its oscillation. Since the radial velocity 
at this point is low, it's likely to remain outside 1/lOth of the virial radius for a significant 
period. The bottom line is that the range of initial kicks that would take the BH to the low- 
density outskirts to cause significant underfeeding, while still allowing it to return quickly 
enough to be "retained" is simply negligibly small. 

Blecha & Loeb (2008) recently performed a more detailed analysis of the orbits of recoil- 
ing SMBHs that include a multi-component halo mass distribution and three-dimensional 
orbits. They report that recoil velocities of between 100 km s~^ and the escape velocity lead 
to significant suppression of the accretion rate, with SMBHs accreting only ~ 10% of its 
initial mass over 10^ — 10^ years of wandering through the halo. We note that (1) in our 
simulations we ignore the longest-wandering BHs through our prescribed retention thresh- 
old; (2) typical kick magnitudes for SMBHs are lower than 100 km s~^ in our simulations, 
a point we explain below; and (3) their prescription of the baryon distribution results in 
a lower central density than our models, as in that paper they are concerned with typical 
galaxies at low redshift, and not with minihalos and protogalaxies. We conclude that pro- 
longed periods of wandering and underfeeding are unlikely to have a significant effect on the 
mass function of high-2; SMBHs as a whole. Oscillations may play a more prominent role 
in the growth history of SMBHs if large kicks are more common (and the retained holes 
tend to be marginally retained), if the halo gravitational potential is more shallow, or if the 
effect of dynamical friction on BH orbits is less than what we have considered in this paper. 
Also, halo triaxiality (Blecha & Loeb 2008) and/or a clumpy mass distribution (Guedes et 
al. 2008) could increase the time that kicked holes take to return to the halo center (or they 
may never return). In principle, this may increase the impact of these oscillations. However, 
in practice, the inner, baryon-dominated regions of the galaxies are likely close to smooth 
and spherical. The BHs that are not ejected in our models typically do not leave these 
regions and so will not be subject to large changes in their orbits from these effects. 

The results presented thus far seem to paint a relatively simple set of relevant parameters 
for SMBH assembly. There is the seeding function /seed, which governs the BH merger rate 
and therefore to a large extent the LISA event rate. The event rate also depends on the 
time-averaged accretion rate, parameterized here by the duty cycle /duty, and the initial 
seed mass Mgeed, as these set the evolution of the observable mass spectrum. The seeding 
prescription and /duty determine the mass function almost entirely if /seed <^ 1- If /seed ^ 0.1, 
then the spins of the BH binary play a role in setting the recoil speeds and the subsequent 
evolution. 

Once typical spin and mass parameters of merging SMBHs are determined by LISA, 
either by direct measurement, or perhaps by extrapolating from detections at lower redshift. 
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combined with the event rate and what is known about the upper end of the mass function, 
this information is hkely sufficient to give at least a strong indication on the typical SMBH 
growth rate and initial seed mass. 

Our simulations above also confirm a result reported by Volonteri & Rees (2006), namely 
that SMBHs form primarily through repeated mergers of the most massive SMBHs merging 
with less-massive (SM)BHs. This is because the gravitational rocket speeds decrease rapidly 
as the mass ratio q decreases. The first few BHs that "outweigh" their neighbors - be it 
through being endowed with more mass at birth, accreting faster or being fortunate enough 
to survive the first mergers - will be less likely to be ejected from their host halos. This 
survival advantage becomes a runaway effect, as each subsequent merger will result in a 
lower value of q for the next merger. 



3.2. Constraints on the SMBH Mass Function 

Now that we have a first-glance grasp of the assembly parameters and their most basic 
observational characteristics, we can turn to identifying actual candidate models for the 
formation of m > lO^M© SMBHs before z ~ 6. 

The suite of models discussed above has demonstrated that there are several feasible 
ways to build the SMBHs. These models have so far focused only on the number density 
of m > lO^M© SMBHs, and ignored indirect constraints that exist on the mass function 
at lower masses. In particular, the total mass density in SMBHs with masses in the range 
IO^Mq <'m< IO^Mq in the local Universe has been estimated by several authors, who find 
~ 4 X IO^Mq Mpc^^ (to within a factor of ~two; see, e.g., the recent paper by Shankar et 
al. 2009a and references therein). Furthermore, a comparison of the locally observed mass 
density with the mass density inferred from accretion by the evolving quasar population 
suggests that ~ 90% of the total local SMBH mass density is attributable to quasar accretion. 
This implies that the total SMBH mass density increased by a factor of ~ 10 between z ~ 6 
and the local Universe (Yu & Tremaine 2002, Haiman et al. 2004; Shankar et al. 2009a), 
which then places an indirect constraint on the SMBH mass function, down to m ~ IO^Mq, 
at z = 6. 

To be specific, we set the following upper limit on the expected total z ^ 6 SMBH mass 
density in m > IO^Mq SMBHs: 

PsMBH,5+(2; = 6) ~ 0.1 X psMBH,6+(2; ?« 0) ~ 4 X lO^Mo Mpc-^ (22) 

That is, we assume that the total mass density of SMBHs with mass m > IO^Mq at 2; = 6 is 
at most 10 percent of the total inferred mass density of SMBHs with mass m > IO^Mq in the 
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local universe. The major caveat to making such an expectation is that we assume that the 
low-mass end of the BH mass function has grown by a factor of 10, and that high-redshift 
quasar luminosity surveys have sufficiently accounted for selection effects of any SMBHs that 
may be hidden by inactivity or by being too dim. 

The analysis that follows below is similar to that of Bromley et al. (2004), who consid- 
ered the upper limit to the z ^ 6 SMBH mass density in weighing the plausibility of various 
z = 6 quasar BH assembly models. The main difference is that Bromley et al. (2004) did not 
consider the gravitational recoil effect. Adding in this recoil makes the problem significantly 
worse. This is because the recoil necessitates more optimistic assumptions in order to build 
up the ~ IO^Mq SMBHs, which tends to increase, by a large factor, the predicted number of 
lower-mass ~ IO^Mq SMBHs, which arise later, and whose growth is therefore less sensitive 
to the kicks. In other words, the kicks preferentially suppress the abundance of the most 
massive SMBHs which arise from the earliest seeds; as a consequence, we predict steeper 
SMBH mass functions than the models considered in Bromley et al. (2004). We performed 
a series of assembly calculations explicitly without recoil, and find that these indeed pro- 
duce mass functions with a fiatter slope and a higher normalization, owing to greater SMBH 
masses at the high-mass end, and a higher occupation fractions at all masses. If there are no 
kicks, less optimistic assumptions are required to produce the SDSS quasars, and the overall 
mass density is lower in no- kick scenarios that produce the minimum number of > 10^ Mq 
SMBHs. 

The basic result we find is that among the models presented thus far, all of those that 
match the SDSS abundance of the most massive SMBHs overshoot the above value by two 
or more orders of magnitude. One possible solution to this apparent problem is that there 
is no problem at all: we have simply set our constraint too severely. Perhaps not all IO^Mq 
SMBHs ai z ^ 6 have evolved to become m > IO^Mq BHs in galactic nuclei in the local 
universe. Still, there is no obvious way to "hide" these low-mass SMBHs locally, and the 
over-prediction in the SMBH densities in our simulations are very large: we find that in our 
fiducial models, we typically need to reduce the population of SMBHs in the mass range 
IO^^^Mq by a factor of > 100, and those in the lO^'^MQrange by a factor of > 10. 



3.3. Successful Models I: BH Seeds Stop Forming Early 

One logical solution, and one that has been suggested by Bromley et al. (2004), is 
to simply stop making seeds below some cutoff redshift. The earliest seeds contribute the 
majority of the total mass of the most massive SMBHs, and late-forming seeds tend to 
end up in lower-mass SMBHs. In principle, then, by introducing a cutoff redshift for seed 
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formation, one can suppress the low-mass end of the mass function, while still allowing for 
the formation of the most massive SMBHs. Such a model would also be in line with our 
physical understanding of seed BHs. We know Pop III stars stopped forming relatively early 
on in the universe, with the halt in production being due to trace metal contamination (e.g. 
Omukai et al. 2008), radiative feedback from the UV and/or X-ray background during the 
early stages of reionization (e.g. Haiman et al. 2000), the higher turbulence of gas in the 
centers of later halos, or a combination of these factors. 

This proposed solution amounts to keeping /duty a constant while allowing /seed to 
evolve with z; in the simplest case, as a step function dropping to /seed = below some 
redshift. Figure [12] shows the fractional contribution {dM/dzsecd)/M from seeds forming at 
different redshifts to the final mass at z = 6 in three different bins of the z = 6 SMBH mass 
function, for two different models that are marginally able to assemble the SMBHs powering 
J1054-I-1024. The model in the left panel has /seed = 10^^, while the model on the right has 
/seed = 1- Note that contributions to each mass bin are normalized to integrate to unity, but 
the two lower mass bins p(lO^~^M0) and p(lO^~^M0) make up the vast majority of the total 
mass density. The formation epochs of the seeds contributing the majority of the mass of 
the different SMBHs are very distinct in the model with lower /seed, but overlap significantly 
for /seed = 1- What accounts for this qualitative difference in the assembly epochs? There 
are two ways to build IO^Mq holes: accretion onto the earliest-forming holes that undergo 
few mergers and grow mostly in isolation, and the mashed-together product of many seeds 
that may have formed later. If the occupation fraction is high, one expects both populations 
to contribute, and therefore a wide distribution for dM / dz^eeA- If the occupation fraction is 
low, SMBHs with many progenitors contributing become rare, and we find that the most 
massive SMBHs can only be formed from isolated seeds in the earliest minihalos, through 
few mergers, and so the dM/dz curves in these models are sharply peaked. 

Our task is to eliminate ~ 99% of the BH mass in the lower-mass bins, while leaving 
most of the m > IO^Mq holes unaffected. By simply examining the normalized dM/dz curves 
in Figure [T2|, one can see that simply cutting off seed formation at an arbitrary redshift for 
the /seed = 1 niodel is not a viable solution to the overproduction problem, as there is no way 
to eliminate the lower-mass SMBHs without eliminating a significant fraction of the IO^Mq 
holes. Cutting off the seed production can produce successful mass functions for models 
with /seed -C 1, but only if the seed cutoff occurs at very high redshifts, typically z ~ 30. 
Essentially, the solution calls for a very brief and early period of seed BH formation, and 
very rare seed formation there after. An unfortunate generic consequence of this early cutoff 
redshift is that it quickly chokes off the LISA event rates. 

We also simulated models where seed production continues beyond the cutoff redshift - 
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with the same probability /seed as in the minihalos aX z > Zcut ^ in halos with Tvir > lO^K. 
Such halos could continue to form BHs if heating by the UV background is the primary 
mechanism for seed suppression, as they are able to shield their central gas from the UV 
radiation (Dijkstra et al. 2004). We find that models with such a partial cutoff overproduce 
the SMBH mass density if /seed ^ 10"*^. This result implies that either the initial occupation 
fraction is very low (it is still possible to make the SDSS BHs with this low /seed; see Table 2 
below), or else some other feedback beyond UV radiation, such as metal-enrichment, stops 
seed BHs from forming in all halos at z < Zcut, even in the rare, more massive ones. 

In short, the requirement in models in which the duty cycle is a constant, but seeds stop 
forming suddenly below some redshift, can be simply summarized as follows: the only way to 
build SMBH mass functions that satisfy observational constraints and indications at both the 
low-mass and the high-mass ends is from extremely rare seeds that form during a brief and 
very early epoch. We also note that these models are attractive because (i) there are physical 
reasons for the seeds to stop forming below some redshift, and (ii) there is independent 
empirical evidence, from constraints on the reionization history from WMAP measurements 
of the cosmic microwave background polarization anisotropies, that the ionizing luminosity 
in high-redshift minihalos was suppressed by a factor of > 10 (Haiman & Bryan 2006). 

We present in Table 2 the parameters for four such successful models. While it is not 
computationally feasible to search the entire parameter space for such models, we present 
two typical examples for both the Pop-HI-remnant and direct -collapse seed BH scenarios, 
/seed is low in each of these examples, as we have argued above that they must be. Although 
we have listed the spin prescriptions, they are relatively unimportant because seeds are rare 
(see Section 3.1 above). Given a particular value for /seed, the only free parameters are the 
accretion rate /duty and the seed cutoff redshift z^ut- For the Pop HI models, we find in the 
range /seed > 10^"^ that seeds must typically stop forming at z^nt > 30, with a lower cutoff 
Zcut ~ 20 for the direct-collapse models. An important conclusion is that in each of these 
models, GW events are too rare (< 10~^/dz/yr for z < 30) for LISA detections beyond 
2; ~ 6 to occur within the mission's lifetime. 

The mass density of ejected BHs is exceedingly low in each of the successful scenarios, 
< 1O~^M0 Mpc~^. Because the total number of seeds is small, so are the number of ejected 
holes. We only give an upper limit here, because the ejected holes are too rare for us to give 
a robust value given the statistical limitations of our "halo cloning" method for populating 
the entire halo population aX z = 6 (the mass density in ejected BHs can be large in models 
with large /seed! see below). 

These models represent the simplest scenarios for SMBH formation, requiring a very 
brief period of seed formation and a prolonged period of accretion at rates comparable to 
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the Eddington rate, and consequentially they represent the most pessimistic predictions for 
LISA'S observational prospects. 



3.4. Successful Models II: Feedback Adjusted to Maintain m-a Relation 

While the suppression of BH seed formation is an attractive possibility that fits con- 
straints on the z = Q SMBH mass function, it is clearly not unique. One alternative solution 
to the over-production problem is to simply reduce the accretion rate of lower-mass BHs at 
lower redshifts. This is again physically plausible: accretion could be choked off as a result 
of the baryonic gas being churned into stars, being heated and dispersed by reionization 
feedback, or through self-induced negative feedback where the BH's own accretion-powered 
radiation stops the gas supply. Rather than try to model such a time- (and probably mass-) 
dependent mass accretion scenario, we examined several model variants, in which BHs are 
allowed to accrete just enough mass to match the value inferred by the m-a relation between 
BH mass and host halo velocity dispersion. That is, at each timestep t ^ t + At, all BHs 
were assumed to grow in mass hj m ^ m + Am such that the m — M relation is satisfied 
at the new host halo mass and redshift. However, whenever this requires super-Eddington 
growth, i.e. if (m + Am)/m > exp(At/tEdd), then Eddington growth is applied instead. 
The main additional assumption here is that the m-a relation remains valid at all redshifts 
(which is at least consistent with a comparison between the evolution of quasars and early- 
type galaxies at < 2; < 6; Haiman et al. 2007; Shankar et al. 2009b). As above, we adopt 
Equation (!20|) as our extrapolated m/M relation. 

The BHs in these models form as lOOM© seeds, and, given their host halo mass and 
redshift, accrete to match this relation as closely as possible without exceeding the Eddington 
accretion rate. If the simulation completes with the mean BH accretion rate well below the 
Eddington rate at all redshifts, then it is consistent with satisfying the Eddington accretion 
rate and the m — M relation inferred by Equation fl20l) . As we shall find below, in our 
models the maintenance of the m-cr relation does not typically require that the Eddington 
accretion rate is saturated (see Figure [T5] below) as long as the /seed ^ 10~^. If both the seed 
mass and the seeding fraction are low, it is increasingly difficult to satisfy Equation (l20l) at 
higher redshifts while simultaneously satisfying the Eddington upper limit. We find that for 
/seed ^ 10~^, accretion must saturate at the Eddington rate for much of 2; > 15 until the 
extrapolated m,-a relation is satisfied, with the mass function falling below this relation at 
earlier stages of growth. 

Recoil velocities are calculated with the spin magnitudes chosen uniformly between 0.0 
and 0.9. As with our previous models, we run simulations where the spins are either randomly 
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oriented or completely aligned with the angular momentum vector of the binary orbit. This 
class of models in effect represents the most optimistic LISA expectations, as it allows us 
to keep numerous seeds, while simply adjusting the accretion rate, as described above, to 
keep the mass function within bounds. Note that the recoil speeds are also minimized by 
our choice for the spin alignment. 

We show the mass functions and occupation fractions for this model in Figure [131 ^oi 
three different values of the cutoff redshift below which new seeds are not formed, Zcut = 0, 
12 and 18. Note that it is still possible to form the SDSS quasar-SMBHs via Eddington- 
limited accretion hj z ^ 6 even if seeds form in just 0.1% of all 1200K halos and only 
before z = 18. We do not plot the m/M relation, as it is satisfied in the form of Equation 
f l2U]) in all cases shown here, by construction. These models also satisfy, by construction, 
the upper limit for the SMBH mass density. In all of the models shown in Figure [T31 
PSMBH,5+(^ = 6) > 1.3 X lO^Mo Mpc-3 if /seed > 10"^ In the /seed = 10"=^ modcls, we 
find psMBH,5+(-2 = 6) ~ 1.0 X lO^M© Mpc~^. The difference in psmbh is due to varying 
occupation fractions at the low end of the halo mass function. The mass functions in Figure 
[13] have a much shallower slope overall than those in Figures [5] and [HI For the mass functions 
shown earlier, the steeper slopes were due to the ratio m/M increasing with decreasing host 
halo mass; the observed m/M relation has the opposite trend. 

We show the LISA event rates in these alternative models in Figure [TH Most signifi- 
cantly, we note that in the /seed = 1 versions of these successful models, the LISA event rate 
can be as high as 30 yr~^. (Note that this number can be even higher if seeds can form in 
minihalos down to a virial temperature that is significantly lower than our assumed fiducial 
value of 1200K.) The rate is somewhat suppressed when compared to the earlier Pop III seed 
models (Figure [Tj) that exceeded realistic indications on the SMBH mass density, because 
the massive BHs in the LISA band are rarer due to the more modest growth rates. We draw 
the attention of the reader to the apparent independence of the detection rate on the seed 
fraction and seed formation cutoff in the cases where /seed ^ 10"^ and Zcut ^ 12. Because BH 
ejections probabilities are lower in these models when compared to the constant-accretion 
scenarios of Figures [7] and [TO], and because the m/M ratios are the same function of M 
in all models shown, the LISA rates saturate and converge once the occupation fraction 
in the most massive halos approach unity. Because the T = 1200K halos form in greatest 
abundance at z ~ 20, the Zcnt = 12 case hardly differs from the case with no seed cutoff. 

A key characteristic of any SMBH assembly scenario is the balance between growth 
through BH mergers and growth through gas accretion. As discussed above, the two must 
strike a balance such that they are able to account for the most massive observed quasar- 
SMBHs at 2; ~ 6, while also not exceeding the total observed SMBH mass density. In Table 
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3, we illustrate the relative importance of mergers vs. growth in the models presented in this 
paper: the four successful constant-accretion models from Table 2; two of the unsuccessful 
constant-accretion models, also in Table 2, which overproduce the universal SMBH mass 
density; and four of the models that explicitly follow the extrapolated m-a relation via 
Equation (!20|) . The values shown in the table (in log^g) are the sum of the initial masses of 
all the seed BHs that enter our merger trees; the total mass of galactic BHs at 2; = 6 13; and 
the total mass of BHs ejected before z = 6. We also calculate the ratio of the total (galactic 
and ejected) BH mass at z = 6 to the total initial seed BH mass, which gives a simple 
measure of the growth through gas accretion. We see immediately the contrast between the 
two types of successful models: the constant-accretion scenario relies on gas accretion for 
much of the growth, typically several orders of magnitude in the total BH mass; where the 
self-regulating models essentially describe the most heavily merger-driven scenarios possible, 
requiring accretion-driven growth of as little as a factor of a few. 

These models also produce a significant population of ejected BHs. Even though ejection 
rates are lower on the whole than our constant-accretion models (compared to the unsuccess- 
ful constant accretion Model X in Tables 2 and 3), two factors contribute to the ejected BH 
mass being comparable to the galactic BH mass at 2; = 6. First, seed BHs are allowed to be 
very common, especially in contrast to the successful constant-accretion Models A through 
D; this results in a far greater number of total merger events, and a high total number of 
ejections despite the lower ejection probabilities. Second, the surviving BHs do not grow 
nearly as rapidly in these models as in the constant-accretion scenarios, so that the ratio 
between the total galactic BH mass and the total ejected mass can remain low (whereas in 
the constant accretion scenarios, retained holes can rapidly outgrow their escaped counter- 
parts.) In the /seed = 1 models, the ejected holes can outnumber and outweigh their retained 
galactic counterparts with mass densities of ~ 3 x W^Mq Mpc~^ and number densities of 
~ 100 Mpc~^. The mass function of the ejected holes is peaked slightly above the seed mass 
(because holes are most likely to be ejected at the earliest stages of their evolution, when 
their host halos are the least massive). All of the ejected BHs are ^ IO^Mq if spins are 
aligned, but in rare instances, SMBHs as massive as ~ IO^Mq are ejected in our models 
with randomly oriented spins (the ejected SMBHs with masses above m > IO^Mq have a 
very low number density, SMBHs of ~ 4 x 10~^ Mpc~^, even in the model with the most 
ejections (random orientation, no cutoff redshift, /seed = !)• 



'"^Note that the values in Table 3 include BHs of all masses equal to and above the seed mass, where we 
have considered only those with m > IO^Mq in computing the universal "SMBH" mass density in previous 
sections. Also note that we do not generate trees for DM halos with M{z = 6) < 10® Mq, and throughout 
this paper we do not account for any BHs that may reside in such halos. 
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These self-regulating accretion models work by adjusting the BH accretion rates ac- 
cording to the mass growth of their host halos. We plot the accretion rates in units of 
the Eddington rate for this new set of models in Figure [13 We do so for the Zcnt = 12 
case with four combinations for /seed and spin alignment, and for three different BH mass 
ranges: IO^Mq < m < IO^Mq, 10^ Mq < m < 10^ Mq, and m > 10^ Mq. Note that the 
accretion rates must be slightly higher if BH binary spins are randomly oriented, in order 
to compensate for the higher ejection rates. Similarly, accretion rates are higher if seeds 
are less common, in order to compensate for the reduced merger-driven growth. For the 
models shown, the duty cycle (the mean accretion rate in Eddington units) for the most 
massive SMBHs converge to ~ 0.2 at z ^ 6, though it can be as high as > 0.5 at 2; > 8 if 
merger-driven growth is hindered by low occupation fraction or recoil-induced ejections. 

We note that similar merger-tree models tracking SMBH growth have been published. 
For example, Koushiappas et al. (2004) presented a similar model where the SMBHs are 
assembled primarily through mergers of directly-collapsed halo cores. Bromley et al. (2004) 
also presented SMBH assembly model wherein gas accretion activity was triggered by major 
mergers of the BHs' host halos; in their model, a set fraction of the baryonic mass of the host 
was fed to the BH at each major merger. Their prescription (albeit without gravitational 
recoil) successfully produced the most massive SDSS SMBHs before redshift z = 6 without 
overproducing the mass density. In general, this type of assembly model is fairly easily tuned 
to broadly reproduce the m — a relation, as the parallel mass growth of BHs and their host 
halos is built in. 



4. Conclusions 

In this paper, we have attempted to map out plausible ways to assemble the > IO^Mq 
SMBHs that power the bright redshift z ^ Q quasars observed in the SDSS, without over- 
producing the mass density in lower-mass (~ 10^~'^Mq) BHs. We also computed the event 
rates expected for LISA in each of the successful models. 

The physical effects governing SMBH assembly depend on the answers to four basic 
questions: (1) how common are the initial BH seeds; (2) how much mass in gas do they 
accrete, and therefore how much they contribute individually to the final SMBH's mass; (3) 
how often do they merge; and (4) what happens to SMBH binaries when they do merge? 
Currently, we do not have empirical constraints to offer definitive answers to any of these 
questions. However, we are capable of predicting the final outcome, starting with a set 
of assumptions for the underlying physics. Our trees-plus-orbits algorithm simulates the 
formation history of SMBHs and the subsequent detection rate expectations for LISA by 
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isolating and prescribing answers to the above four questions. It is a powerful simulation 
tool, as it can incorporate a detailed modeling of individual physical prescriptions without a 
significant increase in the computational load, as long as the prescriptions can be described by 
fitting formulae, tabulated in a lookup table of reasonable size or summarized in a statistical 
manner. 

Using this tool, we have surveyed a wide range of candidate assembly models, and 
reported on common and distinguishing traits in the resulting SMBH mass functions and 
the corresponding LISA detection rates. In particular, we have shown that SMBHs can form 
in a manner consistent with other observational evidence either through the rapid growth of 
rare, massive seeds, or through ultra-early production of numerous Pop-Ill remnant seeds, 
provided these seeds stop forming below a redshift z^ut ~ 20 — 30. We reach the pessimistic 
conclusion that these scenarios do not produce any detectable LISA events at z > 6 (at 
least not in a few year's operation). An alternative model, in which we assume that the 
extrapolation of the local m — M relation holds at all redshifts (e.g. due to internal feedback), 
on the other hand, can produce up to ~ 30 LISA events per year, with a characteristic mass 
spectrum. 

Our major findings can be summarized more specifically as follows: 

• SMBHs must be continuously surrounded by dense gas that was able to cool at the 
centers of DM halos. Feeding holes with the low-density gas in DM halos whose gas 
was unable to cool does not allow for high enough accretion rates to explain the SDSS 
quasar BHs. 

• If embedded in dense gas nearly continuously, ~ IOOMq seed BHs can grow into the 
SDSS quasar BHs without super-Eddington accretion, but only if they form in miniha- 
los at 2; > 30 and subsequently accrete > 60% of the time. However, these optimistic 
assumptions, required to explain the SDSS quasar BHs, overproduce the mass density 
in lower-mass (fewxlO^Mo < Mbh < fewxlO'^Mo) BHs by a factor of 10^ - 10^ We 
find that two conditions need to be satisfied to alleviate this overprediction: the ini- 
tial occupation fraction of seed BHs has to be low (/occ ^ 10^^), and new seeds must 
stop forming, or the seeds must accrete at severely diminished rates or duty cycles, 
at z ^ 20 — 30. We argued that models in which BH seeds stop forming at z ~ 20 
are attractive because there are physical reasons for the seeds to stop forming below 
some redshift (such as metal pollution and/or radiative feedback that suppresses pop- 
III star formation), and because there is independent empirical evidence, from WMAP 
constraints on the reionization history, that star and/or BH formation in high-redshift 
minihalos was suppressed by a factor of > 10 (Haiman & Bryan 2006). 
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• The simplest SMBH assembly scenarios, which have constant accretion rates, but in 
which BH seed formation stops abruptly at some redshift, and which meet constraints 
at both the high-mass and low-mass end of the z = 6 SMBH mass function, predict 
negligibly low LISA event rates. The reason for the low rates is as follows: in these 
models, the BHs that grow into the most massive, highest-redshift quasar-SMBHs 
accrete at the same (exponential) rate as all the other BHs, typically resulting in a vast 
overproduction of massive {m ~ 10^ Mq) holes. In order to offset this overproduction, 
seeds must be made very rare, and this diminishes the LISA rates. It is difficult to 
envision a scenario for high (> 10 per year per unit redshift) detection rates unless 
a vast number of SMBHs in the lO^'^M© range lurk in the universe at all redshifts, 
which the current electromagnetic surveys have missed. 

• A different class of successful models, in which the SMBH masses are self-regulated by 
internal feedback, can evade this constraint, and produce LISA rates as high as 30 yr~^. 
The key difference in these models that predict higher LISA rates is that the SMBH 
growth is driven by a large number of seed BHs and far lower gas accretion rates than 
those required in the constant-accretion models. The majority of these events occur 
aX z ^ 6 and in the low end (10^ — 10'^ Mq) of LISA^s mass range for detection. Also, 
for these models we find the ejected BH mass density can exceed that of the galactic 
BH population aX z = 6. Most ejected holes are expected to have masses similar to the 
seed mass, but an ejected BH can be as massive as ~ IO^Mq if large recoil velocities 
are allowed (e.g. if spins are not always aligned with the orbital angular momentum 
of the binary) . 

In addition to the above, our modeling reveals a number of interesting aspects of SMBH 
assembly. We find that in the successful models the initial seeds are rare, and the most 
massive SMBHs grow primarily from the few 'lucky' early seeds that avoided ejection due to 
kicks. The precise assumptions regarding the kick velocity distribution (such as the assumed 
spin orientations or the resulting oscillation of the BH) tend to have only a modest effect on 
the final results in these models. This is because, at least in our simple prescription, BHs 
either return quickly to the gas-rich nucleus or are left wandering in the outer regions. 

Our results suggest that LISA will be capable of narrowing the field of plausible SMBH 
assembly models from the raw event rate, even without detailed measurements of the binary 
spins or mass ratios. The spin and mass ratio measurements will further constrain the 
evolution of SMBH properties. While the component prescriptions explored in this paper 
are admittedly crude, exercises similar to the one performed in our study will be crucial in 
understanding the limits and possibilities offered by LISA, and ultimately to interpret the 
detected LISA events. The scarcity of empirical constraints on the various pieces of physics 
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that determines the SMBH growth leaves us with a large range of "plausible" scenarios and 
free parameters for SMBH assembly. 
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Table 1: Masses and Quantities of Simulated DM Halos. 



M,,<M< Mhi 



hi 



logioA^bin logio(Mbin) A^si 



W^ 



bin 



8.0 < logioM< 8.5 

8.5 < logio M < 9-0 

9.0 < logioM< 9.5 

9.5<logioM< 10.0 

10.0<logioM< 10.5 

10.5<logioM< 11.0 

11.0<logioM< 11.5 

11.5<logioM< 12.0 

12.00 < logio M < 12.5 

12.50 < logio M < 12.85 

12.85 < logio M < oo 



12.33 


8.22 


50000 


4.28 X 10^ 


11.78 


8.72 


27000 


2.23 X 10^ 


11.20 


9.22 


15000 


1.06 X 10^ 


10.57 


9.71 


9000 


4.13 X 10^ 


9.87 


10.2 


5000 


1.48 X 10^ 


9.08 


10.7 


2700 


4.45 X 10^ 


8.14 


11.2 


1500 


9.20 X 10^ 


6.98 


11.7 


900 


1.06 X 10^ 


5.50 


12.1 


500 


632 


3.48 


12.6 


303 


10. 


1.60 


12.9 


40 


1.0 



For each mass bin, BH assembly was simulated by creating a merger tree for A^sim 
Monte-Carlo halos, and the results were multiplied by the weighting factors Wun to repre- 
sent the Press-Schecter mass function for DM halos at z = 6. A^bin = /m'" dN/dM dM is the 
expected Press-Schechter number of halos in each bin, and (Mbin) is the number-weighted 
mean halo mass in each bin. 



Table 2: Properties of four successful (A-D) and two failed (X and Y) models for SMBH 
growth 



Model 


^^seed 


-^ seed 


Jsccd 


/duty 


spin 


^cut 


PSMBH,5+(2; = 6) 


A 


2OOM0 


1200K 


10-4 


1 


aligned 


25 


3.4 X 10*Mq Mpc-3 


B 


IOOM0 


1200K 


10-2 


0.95 


aligned 


28 


5.1 X 10^ Mo Mpc-3 


C 


lO^Mo 


1.5 X lO^K 


10-4 


0.6 


random 


13 


6.2 X 10*Mq Mpc-3 


D 


2 X IO^Mq 


1.5 X lO^K 


10-2 


0.55 


aligned 


18 


7.0 X 10*Mq Mpc-3 


X 


IOOMq 


1200K 


1 


0.8 


random 





2.9 X IO^Mq Mpc-3 


Y 


lO^Mo 


1.5 X lO^K 


10-3 


0.6 


aligned 





1.1 X 10'^ Mq Mpc-3 



The above shows parameters for four models that (1) have constant accretion rates of 
/duty times the Eddington rate; (2) produce by 2; = 6 SMBHs massive enough to power the 
SDSS quasars; and (3) do not overproduce the overall SMBH population. In Models A and 
B the seed BHs are Pop III remnants, and in Models C and D the seeds are formed through 
direct-collapse in more massive halos. Models X and Y are unsuccessful models that barely 
produce the m > IO^Mq SMBHs by 2; = 6 but also far overproduce the lower-mass SMBHs. 
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Table 3: The total seed, total ejected, and total retained BH masses 



Model 






^^seed, tot 


mGN(^ = 


6) 


^^cjected, tot 


"^BH, tot/'"^seed, tot 


A 






9.3 


16.0 




6.2 


6.7 


B 






10.1 


16.1 




7.7 


6.0 


C 






12.7 


16.2 




6.4 


3.4 


D 






13.2 


16.3 




6.2 


3.1 


X 






16.0 


20.0 




19.1 


4.0 


Y 






14.6 


17.5 




12.7 


2.9 


m-a, Zcut = 0, /seed - 


= 1, 


aligned 


15.7 


15.7 




16.0 


0.48 


m-a, Zcnt = 0, /seed = 


= 1, 


random 


15.7 


15.7 




16.3 


0.70 


m-a, Zcut = 0, /seed = 


10- 


^, aligned 


13.7 


15.6 




12.9 


1.9 


m-a, Zcnt = 18, /seed 


= 1 


, aligned 


14.4 


15.7 




14.7 


1.34 



The decomposition of the final SMBH mass into the contributions from the initial stel- 
lar seed BHs and subsequent gas accretion. Masses are in log^g; ^^^ i^ units of Mq. The 
columns, starting from the second and from left to right, show: the total initial seed mass; 
the total SMBH mass retained in nuclei aX z = 6 ; the total ejected BH mass, and log^g of ^^^ 
ratio of the total (i.e. the sum of the ejected and nuclear) z = 6 SMBH mass to the initial 
seed mass. The last ratio is a measure of the total growth due to gas accretion. The first 
four rows (A-D) show values for the four successful models. Also shown for comparison are 
values from two of the unsuccessful, constant-accretion models (X and Y) that overproduce 
the total BH mass function. The model parameters for Models A, B, C, D, X and Y can be 
found in Table 2. The models where the m-a relation is enforced by hand can grow primarily 
through mergers, with gas accretion adding as little as a factor of a few to the total SMBH 
mass aX z = 6. If /seed is sufficiently high, they also eject a total mass in low-mass BHs that 
is comparable to the retained nuclear population (most of the ejected holes have a mass near 
the seed mass). 
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Fig. 1. — The Monte-Carlo-generated mass function of progenitors for a IO^^Mq, zq = 6 
parent halo at 2; = 8, 13, 21 and 34. The histogram is the mean number of 100 realizations, 
and the error bars demarcate the Poisson errors. The solid curve is the EPS prediction. 
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Fig. 2. — The kick velocity distribution as a function of the BH binary mass ratio, after 10^ 
reahzations of Equations ([3]-8) at each value of q. The left panel shows kicks for random 
spin orientations, while the right panel shows kicks when both BH spins are aligned with the 
orbital angular momentum. In each case, the spin magnitudes are chosen from a uniform 
random distribution in the range 0.0 < ai^2 < 0.9. The solid curves show the mean, the 
dashed curves show the l-a range, and the dotted curve gives the maximum value generated 
in the 10^ realizations. 
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Fig. 3. — Examples of the radial motion of a recoiling black hole, for three different mass 
profiles for the host halo. The black curve shows the motion in a pure NFW halo; the halo 
of the blue (dotted) curve assumes that the host galaxy has a DM halo with an NFW profile 
and a corey gas component; the red (dashed) curve assumes a DM halo and a cuspy r~^'^ 
power law gas profile. In all cases, the halo mass is IO^Mq, the BH mass is IO^Mq, the 
redshift is z = 20 and the kick velocity is 100 km s~^. 
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Fig. 4. — The maximuin possible accreted mass by redshift z = 6 for a seed BH born with 
a seed mass mgeed in a halo with virial temperature T^i^ = 1200K at redshift z. If the BH 
is always surrounded by a steep gas profile, with a power-law cusp, the growth remains 
Eddington throughout (which appears as a straight line in this log- log plot). If the gas 
profile has a flat core (as in the TIS profile), the central density is initially insufficient to 
feed the BH at the Eddington rate, resulting in much slower growth. The cuspy and corey 
gas distributions are demarcated by thick and thin lines, respectively, and different seed BH 
masses are shown in different colors (line styles). 
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Fig. 5. — The comoving number densities of SMBHs in different mass bins at redshift 
z = 6. Colored figures are available in the online version. The 24 different models shown 
in the figure assume different parameter combinations as follows. The columns, from left 
to right, adopt /seed = 10~^, 10~^, 10~^, 1. The top row is for simulations with random 
binary spin orientation, and the bottom row is for spins aligned with the binary's orbital 
angular momentum. Time-averaged accretion rates are distinguished by color: black (solid, 
/duty = 1), blue (dot, /duty = 0.8), and green (dash-dot, /duty = 0.6). The numbers in the 
upper-right corners represent \ogiQ[p,/{MQ Mpc~'^)] for each model, in descending order of 
/duty The red (dashed) line demarcates the rough indication for the minimum number of 
z ^ 6 SMBHs in the observable universe with m > 10^'^Mq given the area surveyed by 
SDSS. 
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Fig. 6. — Properties of the SMBH population at ^ = 6 as a function of the halo mass M: 
the percentage of DM halos hosting a central BH (assumed at most to be one BH per halo; 
top rows in both the upper and lower panels) with m > 10~^M , and the mean BH-to-halo 
mass ratio m/M for the halos that do host a BH (bottom rows). Color (line-style) and panel 
schemes are the same as in Figure [H The red (dashed) line is the empirical m/M relation 
extrapolated to z = 6 (see Equation [2U1 Wyithe & Loeb 2003, Ferrarese 2002). Our m/M 
relation has the opposite trend with respect to halo mass from the trend observed in the 
local universe. Note that in some cases, the central BH consumes most of the baryonic mass 
in the halo. 
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Fig. 7. — Expected number of LISA detections per redshift per year due to SMBH mergers 
with binary mass 10'^Mq < {mi + m2)(l + z) < 10"^ Mq. The color, line-type and panel 
schemes are the same as in Figures [5] and [61 
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Fig. 8. — The z = 6 SMBH mass function in the direct collapse scenarios, with mgeed = 
IO^Mq and Tseed = 1-5 x lO^K. Color and panel organization for accretion rate, seed fraction 
and spin alignment is the same as in Figure [51 
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Fig. 9. — The SMBH occupation fraction and the m/M ratios in the direct collapse models. 
Refer to Figure [6] for color and panel organization. The dotted line is the extrapolated 
empirical m/M relation. 
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Fig. 10. — LISA event rates in the direct collapse scenarios. Refer to previous figures for 
color and panel organization. Note the significant reduction in the event rates relative to 
the pop-III seed models, owing to the smaller number of seed-forming halos in the tree. 
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Fig. 11. — Properties of the SMBH population under several variants of our fiducial models. 
All models plotted have /duty = 0.65, /seed = 1, and aligned binary spins, and modify a single 
aspect of the basic fiducial model prescription, as labeled: the gas density is an isothermal 
power-law instead of the fiducial p oc r~^-^ (dark blue, dotted curve); spin magnitudes are 
near-maximal at ai2 = 0.9 (green, dot- short- dash curve); halos are allowed to merge and 
form BH binaries irrespective of their mass ratio (yellow, dashed curve); Pop 111 stars do not 
blow out the gas from their host halos prior to leaving a seed BH (magenta, dot-long-dash 
curve); and recoiling BHs wandering in low-density regions continue to accrete efficiently 
(light blue, long-short-dash curve). We also ran simulations with a corey, TIS gas profile for 
the host halos. We were unable to produce SMBHs of even > lO^M© in those models even 
when prescribing the most optimistic values for the other assembly parameters (therefore 
results from these models are not shown in the figure). 
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Fig. 12. — The fractional contribution to the mass of 2; = 6 SMBHs from IOOMq seed BHs 
that form at different redshifts Zgeed- The contributions are computed in three mass bins 
oi z = 6 SMBHs: 10^ Mq < m < lO^M© (black, solid lines), WMq < m < 10^ Mq (red, 
dotted), and m > 10^ Mq (magenta, dash-dot). The most massive z = 6 SMBHs arise mainly 
from the earliest 1200K progenitors of the most massive halos {z > 20), whereas the seeds 
contributing most of the mass of lower-mass SMBHs formed later {z < 20). 
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Fig. 13. — The SMBH mass functions and occupation fractions at 2; = 6 for various models 
satisfying the m-a relation. These mass functions are much less steep in comparison to those 
shown in Figures [5] and [HI Occupation fractions are higher than those in Figures [6] and [H as 
the BH binary mass ratios are generally lower compared to those models, resulting in less 
powerful recoil kicks. 
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Fig. 14. — The LISA rates in models that satisfy the m-a relation at all redshifts as closely 
as possible without exceeding the Eddington accretion rate. Note that the rate saturates at 
30 yr~^ at 2; > 6 for high seed fractions and low redshifts for seed cutoff. 
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Fig. 15. — The accretion rate (thick lines) in Eddington limits in models that match the 
m-a relation (Equation [20!) at all redshifts, for different ranges of BH masses. We also 
plot the merger rates (thin lines) in units of mergers per halo for different halo mass bins 
corresponding to Mi + M2 = W^m. Note that masses are defined instantaneously at each 
redshift interval, rather than tracking the histories of the z = 6 holes and halos. 



